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FOREWORD 


A  mission  of  the  Defense  Mapping  Agency  (DMA)  requires  determining  geodetic-quality 
absolute  positions  autonomously  and  in  areas  far  removed  from  other  sites.  The  use  of  Differential 
Global  Positioning  Systems  (DGPS)  techniques  or  relative  positioning  is  sometimes  precluded. 
Operating  in  an  unclassified  mode,  keyed  receivers  can  output  Precise  Positioning  Service  (PPS) 
navigation  solutions.  Because  they  are  computed  in  real  time,  PPS  solutions  necessarily  use  the 
broadcast  ephemerides  and  consequently  are  subject  to  ephemerides  and  satellite-clock  prediction 
errors.  The  PPS  solutions  can  be  recorded,  but  the  broadcast  ephemerides  and  the  observations, 
whai  corrected  for  Selective  Availability  (SA),  cannot.  This  report  demonstrates  a  method  whereby 
the  field-recorded  PPS  solutions  can  be  reprocessed  with  the  more  accurate  precise  ephemerides. 
Reprocessing  with  the  precise  ephemerides  offers  improved  quality  in  navigation  solutions,  but  does 
not  require  that  the  SA-corrected  pseudorange  or  phase  observations  be  saved.  This  method  is  called 
Precise  Absolute  Navigation  (PAN). 
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MOTIVATION 


The  purpose  of  the  Precise  Absolute  Navigation  (PAN)  development  was  to  investigate  the 
feasibility  of  improving,  by  postprocessing,  the  accuracy  of  a  set  of  typical  Precise  Positioning  Service 
(PPS)  navigation  solutions.  Real-time,  stand-alone  navigation  solutions  necessarily  use  the  broadcast 
ephemerides  as  the  source  of  the  satellite  position  and  satellite  clock  estimates,  but  in  certain  cases, 
the  particular  mission  could  benefit  from  more  accurate  navigation  solutions.  If  the  mission  is  not 
time-critical,  delaying  the  use  of  the  navigation  solutions  may  be  acceptable  until  the  postfit  precise 
ephemerides  are  available.  In  this  case,  significantly  better  results  may  be  expected. 

With  a  precise  ephemerides,  static  absolute  position  solutions  can  achieve  sub-meter 
repeatability.^’^  Meter-level  navigation  solutions  have  been  demonstrated  by  postprocessing  Lj 
observations  with  postfit  precise  ephemerides,  clock  estimates,  and  an  ionospheric  model.^ 
Applications  where  precise  absolute  na'vigations  would  be  of  interest  include  cases  where  differential 
or  relative  data  are  not  available,  or  where  the  track  of  a  moving  vehicle  or  ship  is  desired  at  remote 
locations  far  from  other  sites. 

In  the  opCTational  environment,  the  PPS  navigation  solutions  will  be  obtained  from  a  user  in 
the  field.  Presumably,  these  solutions  would  be  obtained  with  a  receiver  that  can  correct  for  the 
effects  of  Selective  Availability  (SA)  and  AntiSpoofing  (AS).  AS  encrypts  and  SA  corrupts  the 
Global  Positioning  System  (GPS)  signals,  making  the  precise  signal  and  the  full  accuracy  of  the 
broadcast  ephemerides  accessible  to  authorized  users  only.  These  receivers  can  display  the  navigation 
solution  in  real  time,  but  the  field-corrected  satellite  ephemerides  and  the  SA-corrected  pseudorange 
and  phase  observations  are  classified  and  therefore  are  not  available  to  the  user.  The  PPS  position 
solutions  are  unclassified  and  can  be  recorded  in  the  field  along  with  the  GPS  time  and  the  satellites 
contributing  to  the  solution.  Even  though  the  original  observations  are  lost,  this  information  is 
enough  to  allow  the  PPS  solutions  to  be  improved  at  a  later  time  with  the  postfit  precise  ephemerides. 


METHOD 


For  this  description  of  the  method,  it  will  be  assumed  that  the  GPS  time  and  the  PPS  solutions 
to  be  improved  are  given,  as  are  the  receiver  type  and  the  details  of  the  navigation  algorithm  it  used. 
The  algorithm  information  may  be  needed  to  establish  which  satellites,  of  those  in  view,  produced  the 
particular  navigation  solution.  Alternately,  the  field  receiver  may  be  able  to  supply  the  satellite 
pseudo  random  noise  (PRN)  numbers  that  were  tracked.  To  compute  the  difference  between  the  two 
ephemerides,  an  independent  source  of  the  SA-corrected  broadcast  ephemerides  and  the  Defense 
Mapping  Agency  (DMA)  precise  ephemerides  and  clock  estimates  are  required.  Finally,  for  testing 
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purposes,  two-frequency,  SA-corrected,  pseudorange  and  phase  observations  from  a  typical  receiver 
need  to  be  added  to  the  requirements. 


STANDARD  POSITIONING  SERVICE  (SPS)  SOLUTIONS 

The  form  of  the  observation  equations  that  produces  the  SPS  navigation  solution  in  the  field 
receiver  can  be  reconstructed  as  in  Equation  (1).  The  original  observation  vector  at  time  4  is  denoted 
by  O  The  ranges  computed  from  the  satellite  ephemerides  and  the  estimate  of  the  current  state 
are  placed  in  the  vector.  The  inherently  nonlinear  problem  is  then  linearized  by  writing  the 
observation  vector  equal  to  the  first  two  terms  of  the  Taylor  series  expansion  for  C^. 

O,-  e,,  (1) 

dX^ 

The  last  term,  represents  the  contributions  from  the  terms  that  are  ignored,  plus  noise.  The 
solution  to  Equation  (1)  is  AX^^,  which  is  found  by  least  squares  in  the  usual  way,  as  indicated  in 
Equation  (2).  The  Wis  the  observation  weight  matrix  that  may  be  incorporated  into  the  formulation. 

AX^  =  (A^WA)-U^W(O^  -  CJ  (2) 


The  matrixyd  in  Equation  (2)  represents  the  matrix  of  partials  written  as  — —  in  Equation  (1).  From 

this  result,  the  SPS  state  estimate  can  be  found  by  adding  the  state  update  to  the  current  state 
estimate;  X^  =  worth  noting  that  if  the  precise  ephemerides  were  available  in  the 

field  along  with  the  broadcast  ephemerides,  the  same  observation  vector,  could  have  been  used 
with  either  of  these  ephemerides  to  produce  a  position  solution. 


PPS  SOLUTIONS 

The  PPS  user  forms  and  solves  the  same  sort  of  equations  as  the  SPS  user.  The  difference 
is  that  corrections  for  SA  are  included.  AS  will  not  be  discussed  since  it  has  no  direct  effect  on  the 
PAN  algorithm.  AS  may  deny  certain  classes  of  receivers  access  to  two-frequency  observations  and 
the  ionospheric  correction.  It  also  denies  the  use  of  the  precision  available  from  P-code  observations. 
The  PPS  user  is  assumed  able  to  recover  the  two-frequency  P-code  and  use  it  in  the  observation 
vector. 


The  PPS  receiver  applies  the  dither  corrections  to  the  pseudorange  and  phase  observations. 
When  applied,  the  original  observation  vector,  is  changed  to  The  receiver  also  applies  the 
epsilon  corrections  to  the  broadcast  ephemerides  messages.  The  corrections  appear  in  the  satellite 
position  part  of  the  vector  and  change  it  to  the  vector.  Note  that  both  of  these  SA-corrected 

vectors  are  required  for  a  field  PPS  solution,  but  neither  can  be  saved  for  later  use.  From  there,  the 
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route  to  the  PPS  solutions  proceeds  in  the  same  fashion  as  for  the  SPS.  The  PPS  observation 
equation  is  the  same  as  Equation  (1)  except  that  the  corrected  vectors  are  used,  as  shown  in 
Equation  (3). 


dC 


hk 


dX, 


(3) 


When  Equation  (3)  is  solved  by  least  squares  in  the  same  manner  as  before,  the  PPS  state  vector 
is  found  from  the  sum:  X^  -  *  ^^tk- 


PPS  WITH  A  PRECISE  EPHEMERIDES 


Suppose  that  precise  ephemerides  were  available  in  the  field  when  the  PPS  solution  was 
computed.  The  observation  equation  would  be  as  shown  in  Equation  (4)  where  the  subscript  p  on 
denotes  that  the  satellite  positions  originate  from  the  precise  ephemerides.  There  is  no  need  for  an 
epsilon  correction  to  the  postfit  precise  ephemerides,  but  the  observation  vector,  on  the  left,  is  the 
same  as  the  dither-corrected  vector  used  in  Equation  (3). 


=  c. 


dX, 


(4) 


Since  this  observation  vector  is  common  to  Equations  (3)  and  (4),  they  can  be  equated.  This  is 
shown  in  Equation  (5). 


dC 


hk 


ax. 


ax, 


(5) 


From  Equation  (3)  the  PPS  field  solution  X^  is  known.  Therefore  it  can  be  used  in  Equation  (5) 
in  place  of  the  initial  estimate  X^.  New  computed  vectors  are  needed  to  reflect  the  change  in  the 
initial  state  estimate;  these  are  indicated  by  the  PPS  superscript.  This  substitution  does  not  change 
the  equality  as  long  as  the  observation  vector  remains  unchanged  and  the  linearity  constraint 
continues  to  be  satisfied.  The  new  form  for  Equation  (5),  incorporating  the  PPS  position,  is  shown 
in  Equation  (6). 


PPS 


dC 


PPS 
hk 


'  bk 


ax. 


PPS 


pt 


dC"‘^ 

— ^AX 

ax,  ' 


+  e 


^PPS 

p^ 


(6) 


If  X^  were  used  instead  of  X^  in  Equation  (3),  and  the  system  solved,  the  new  state  solution  would 
be  near  zero.  This  result  is  usu^y  called  the  adjusted  0  minus  C  and  reflects  the  level  of  observation 
noise  and  residual  unmodeled  effects.  Therefore  the  second  term  from  the  left  in  Equation  (6)  can 
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be  neglected,  because  the  state  vector  is  approximately  zero.  When  this  is  the  case, 

Equation  (6)  reduces  to  a  simpler  form  shown  in  Equation  (7). 


A 

'^bk 


c  r 


dCj^^ 

— 

dx^  ' 


+  e. 


PPS 


(7) 


Equation  (7)  can  be  solved  by  least  squares  for  AX^  in  an  expression  like  Equation  (2). 


AX^  = 


(8) 


Then  the  improved  state,  with  the  quality  of  the  postfit  precise  ephemerides,  is  X^  and  is  found  from 
the  sum:  X^  =  xj^^  +  AX^.  This  is  the  PAN  solution.  The  foregoing  argument  is  equivalent  to 

noting  in  Equation  (3)  that  when  AX^^  is  zero,  the  vector  plus  some  noise  is  equal  to  the 
observation  vector.  Therefore  it  is  possible  to  approximately  recreate  the  lost  observations  from 
knowledge  of  the  satellite  ephemerides  and  the  solution  vector  X^.* 


PAN  EVALUATION 


PAN  assumes  that  PPS  field  solutions  exist,  so  the  first  task  on  the  way  to  verification  of  the 
method  was  to  generate  simulated  PPS  solutions  from  GPS  observations  already  on  hand.  This 
required  that  a  navigation  algorithm  be  developed  and  tested.  The  formulation  for  this  algorithm  is 
attached  as  Appendix  A.  The  solutions  obtained  from  the  navigation  algorithm  were  intended  to 
substitute  for  the  PPS  field  solutions  that  are  the  starting  point  for  the  PAN  algorithm.  The 
navigation  solutions  require  GPS  observations,  and  two  sources  of  existing  data  were  available  as 
described  in  the  following  paragraphs. 

DMA  provided  one  data  set  for  one  week  of  observations  at  static  site  85401  in  March  1995. 
An  AshtechZ12  receiver  was  used.  The  data  interval  was  30  sec  and  both  SA-uncorrected  and  SA- 
corrected  observations  were  provided.  For  the  PAN  application,  the  SA-corrected  data  were  used 
in  the  navigation  algorithm  to  simulate  PPS  solutions  from  a  field  site.  Coordinates  for  this  site  were 
provided  with  the  data,  but  absolute  position  solutions  were  also  performed  on  three  of  the  days  as 
a  check.  The  site  documentation  is  included  as  Appendix  B. 

The  second  set  of  observations  were  selected  from  the  Drive-By  experiment  performed  at  the 
Naval  Surface  Warfare  Center,  Dahlgren  Division  (NSWCDD),  Dahlgren,  Virginia,  in  February  1994. 
Four  Trimble  4000  SSE  receivers  collected  1-sec  data  simultaneously  for  tests  of  Remondi's  on-the- 
fly  (OTF)  algorithm.^  Two  of  these  receivers  were  at  fixed  sites  and  two  traveled  several  laps  around 
a  circuit  on  the  naval  base  with  antennas  attached  to  the  roofs  of  two  vehicles.  The  data  sets  were 
given  to  DMA  in  May  1995  so  that  the  SA  effects  could  be  removed  in  the  DMA  Data  Correction 
Facility  (DCF).  The  SA-corrected  data  are  classified,  and  were  processed  by  the  navigation  algorithm 
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in  a  classified  facility.  The  resulting  unclassified,  simulated  PPS  solutions  were  compared  with  the 
OTF  results,  which  represented  the  true  vehicle  track  with  an  estimated  error  of  a  few  centimeters. 


STATIC  SITE  RESULTS 

The  true  absolute  position  for  85401  was  obtained  relative  to  a  primary  control  station, 
USAF3  1988.  The  errors  in  the  relative  positioning  are  given  as  ±0.02  m  in  each  component,  while 
the  absolute  error  in  USAF3  is  given  as  ±1  m  in  each  component  with  respect  to  the  World  Geodetic 
System  (WGS84)  datum.  These  coordinates  are  used  throughout  this  section  as  the  true  position  for 
85401. 


The  SA-corrected  observations  and  broadcast  ephemerides  were  processed  with  the 
navigation  algorithm  to  obtain  the  simulated  PPS  field  solutions.  The  PPS  solutions  were  then  used 
with  the  DMA  precise  ephemerides  and  clock  estimates  to  perform  the  PAN  solutions  at  each  PPS 
epoch.  The  daily  PAN  averages  showed  a  consistent  1.5-  to  2.0-m  offset  in  the  north  and  -1.0  m 
east,  compared  with  the  truth  (PAN  -  Truth).  Since  this  disagreement  may  be  have  been  due  to  small 
errors  in  the  truth,  three  days’  data  were  processed  with  the  DMA  absolute  positioning  program 
GASP.®  The  results  fi’om  GASP  using  the  SA-corrected  broadcast  ephemerides  are  shown  in  Table  1; 
the  results  using  the  precise  ephemerides  are  shown  in  Table  2.  As  expected,  the  consistency  between 
the  three  days  is  better  with  the  precise  ephemerides.  This  is  quantified  by  the  smaller  standard 
deviations  for  the  latitude,  longitude,  and  height  in  Table  2. 


TABLE  1 .  GASP  SOLUTIONS  FOR  85401  WITH  SA-CORRECTED  BROADCAST  EPHEMERIDES 


DATE 

1995 

LAT 

DELTA 

(m) 

LONG 

DELTA 

(m) 

ELHGT 

DELTA 

(m) 

RSS 

(O)* 

(m) 

X 

DELTA 

(m) 

Y 

DELTA 

(m) 

Z 

DELTA 

(m) 

RSS 

(C)*» 

(m) 

RMS 

(cm) 

NO. 

OF 

OBSV 

% 

REJEC. 

TED 

072 

0.081 

-0.131 

1.306 

1.315 

-0.135 

-0.969 

0.879 

1.315 

4.900 

8060 

1,0 

073 

-0.867 

-2.065 

1.686 

2.803 

-2.072 

-1.850 

0.376 

2.803 

9.100 

7799 

3.0 

075 

0.997 

-1.373 

-1.212 

2.085 

-1.367 

1.575 

0.022 

2,085 

4.500 

8097 

1.0 

MEAN; 

0.070 

-1.189 

0.593 

2.068 

-1.191 

-0.415 

0.426 

2.068 

6.167 

STD: 

0.932 

0.980 

1.575 

0.744 

0.981 

1.778 

0.431 

0.744 

STDM: 

0.538 

0.566 

0.909 

0.566 

1.027 

0.249 

TOTAL  NUMBER  OF  DAYS  PROCESSED:  3 

*  Geodetic 
**  Coordinates 


The  mean  solutions  from  Table  2  move  the  station  location  0.3  m  north  and  -1.03  m  east. 
If  both  of  these  changes  were  adopted,  the  east  PAN  position  error  would  be  nearly  zero,  but  a  1 .5-m 
disagreement  would  remain  in  the  north  component.  The  cause  of  this  north  error  is  unresolved  at 
this  writing.  The  average  PPS  and  PAN  navigation  solutions  for  each  day  are  listed  in  Tables  3  and  4. 
Table  3  includes  the  total  number  of  30-sec  epochs  processed  on  each  day. 


5 


NSWCDD/TR-95/196 


TABLE  2.  GASP  SOLUTIONS  FOR  85401  WITH  PRECISE  EPHEMERIDES 


DATE 

1995 

LAT 

DELTA 

(m) 

LONG 

DELTA 

(m) 

EL  HOT 
DELTA 
(m) 

RSS 

(G)» 

(m) 

X 

DELTA 

(m) 

Y 

DELTA 

(m) 

Z 

DELTA 

(m) 

RSS 

(C)*» 

(m) 

RMS 

(cm) 

NO. 

OF 

OBSV 

% 

REJEC¬ 

TED 

072 

0.184 

-1.049 

-0.010 

1.065 

-1.049 

0.127 

0.137 

1.065 

4.200 

8028 

1.0 

073 

0.432 

-1.028 

-0.609 

1.270 

-1.025 

0.749 

-0.042 

1.270 

4.300 

8038 

0.0 

075 

0.284 

-1.025 

-0.709 

1.279 

-1.023 

0.735 

-0.220 

1.279 

4.100 

8015 

1.0 

MEAN: 

0.300 

-1.034 

-0.443 

1.205 

-1.032 

0.537 

0.042 

1.205 

4.200 

STD: 

0.125 

0.013 

0.378 

0.121 

0.015 

0.355 

0.179 

0.121 

STDM: 

0.072 

0.008 

0.218 

0.008 

0.205 

0.103 

TOTAL  NUMBER  OF  DAYS  PROCESSED;  3 

*  Geodetic 
**  Coordinates 


TABLE  3.  DAILY  STATIC  PPS  WEIGHTED  MEANS  (^M)  AND  STANDARD  DEVIATIONS  (oM) 


m 

071 

072 

073 

074 

075 

O" 

16 

0' 

11 

■ 

n 

0 

n 

0 

n 

G 

m 

B 

B 

0 

B 

0 

B 

o 

E* 

-0.55 

2.27 

-0.44 

2.53 

-0.90 

2.89 

1.07 

1.68 

-0.42 

2.21 

-1.87 

8.69 

1.01 

1.81 

N** 

-0.28 

2.99 

0.23 

3.28 

0.52 

5.12 

0.80 

4.29 

0.56 

3.12 

1.45 

8.07 

-0.52 

2.75 

Ut 

1.09 

6.52 

-0.60 

9.35 

1.14 

8.86 

2.37 

5.65 

1.80 

7.88 

-3.38 

17.26 

0.62 

4.95 

Stt 

23: 

93 

2833 

2833 

1652 

2833 

2826 

719 

*  East  Component 
**  North  Component 
t  Vertical  Component 

tt  Number  of  30-sec  observations  processed  each  day 


TABLE  4.  DAILY  STATIC  PAN  WEIGHTED  MEANS  (pM)  AND  STANDARD  DEVIATIONS  (oM) 


■ 

071 

072 

073 

074 

075 

076 

O'; 

n 

■ 

B 

0 

B 

0 

B 

G 

B 

G 

B 

G 

B 

G 

B 

G 

E* 

-0.78 

0.66 

-0.93 

0.27 

-0.92 

1.88 

-0.90 

0.62 

-0.85 

0.66 

-0.66 

2.04 

-1.12 

0.56 

N** 

1.82 

1.27 

1.77 

1.27 

1.87 

1.36 

1.79 

1.23 

1.88 

1.30 

1.90 

4.01 

1.82 

1.34 

ut 

0.09 

1.88 

0.34 

1.94 

0.17 

3.21 

0.70 

1.80 

0.19 

1.85 

-0.41 

8.29 

0.42 

1.75 

*  East  Component 

**  North  Component 
t  Vertical  Component 
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Figure  1  shows  a  plot  of  the  Geometric  Dilution  Of  Precision  (GDOP)  for  the  navigation 
solutions  obtained  for  85401.  The  elevation  cutoff  angle  was  set  at  15  deg  in  order  to  restrict  the  use 
of  low  elevation  satellites  and  improve  the  quality  of  the  solutions.  The  receiver  tracked  all  satellites 
down  to  10  deg,  but  there  were  frequent  jumps  in  the  solutions  when  data  at  the  lower  elevation 
angles  were  included.  The  GDOP  peaks  at  over  600  during  a  short  span  of  the  day.  The  satellite 
elevation  angles  during  this  period  are  plotted  in  Figure  2.  Just  four  satellites  are  in  view  above 
15  deg  for  about  2300  sec  beginning  at  about  142,200  sec. 


GDOP  FOR  85401 

GPS  WEEK  792,  ELEVATION  ANGLE  CUTOFF  15  DEG 


GPS  TIME  OF  WEEK  (sec) 


FIGURE  1.  GDOP  AT  85401  WITH  A  15-DEG  ELEVATION  CUTOFF 


The  adjusted  observed-minus-computed  residuals  are  the  differences  between  the  observations 
and  the  computed  values  after  the  current  state  solution  is  substituted  for  the  predicted  state.  The 
least  squares  solution  minimizes  the  square  of  the  residuals;  therefore,  the  adjusted  residuals  have 
zero  mean.  The  RMS  of  the  pseudorange  residuals  for  the  PPS  and  PAN  solutions  for  site  85401  and 
day  072  are  plotted  in  Figure  3.  The  jumps  in  the  residuals  are  due  to  new  satellites  being  included 
in  the  solution  or  satellites  dropping  out  of  the  solution.  Only  four  satellites  were  included  in  the 
solution  during  periods  where  the  adjusted  residual  is  zero.  Since  the  two  curves  are  similar,  the 
source  of  the  unmodeled  residuals  must  be  common  to  both  solutions. 

Figures  4  to  17  illustrate  the  differences  between  the  navigation  solutions  and  the  true  position 
for  85401  on  each  of  the  seven  days.  The  seven  even-numbered  figures  show  the  PPS  solutions,  and 
the  odd-numbered  figures  show  the  PAN  solutions.  The  improvement  due  to  the  precise  ephemerides 
is  evident  in  the  PAN  solutions.  The  large  variations  that  remain  in  the  solutions  are  primarily 
associated  with  high  GDOP  periods. 
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FIGURE  2.  SATELLITE  ELEVATION  ANGLES  DURING  PERIOD  OF  HIGH  GDOP 


FIGURE  4.  PPS  SOLUTIONS  FOR  STATIC  SITE  85401,  DAY  071 


FIGURE  5.  PAN  SOLUTIONS  FOR  STATIC  SITE  85401,  DAY  071 
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FIGURE  6.  PPS  SOLUTIONS  FOR  STATIC  SITE  85401,  DAY  072 
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NAVIGATION  SOLUTIONS  FOR  85401,  DAY  073 

SIMULATED  PPS  FROM  THE  BROADCAST  EPHEMERIDES 


East:  mean=-0.90  a=2.89 
North:  mean=0.52  a=5.1 2 
Up:  mean=1 .14  0=8.86 


NAVIGATION  SOLUTIONS  FOR  85401,  DAY  073 

PRECISE  ABSOLUTE  NAVIGATION 


FIGURE  9.  PAN  SOLUTIONS  FOR  STATIC  SITE  85401,  DAY  073 
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NAVIGATION  SOLUTIONS  FOR  85401,  DAY  074 

SIMULATED  PPS  FROM  THE  BROADCAST  EPHEMERIDES 


FIGURE  10.  PPS  SOLUTIONS  FOR  STATIC  SITE  85401,  DAY  074 


NAVIGATION  SOLUTIONS  FOR  85401,  DAY  074 

PRECISE  ABSOLUTE  NAVIGATION 


GPS  TIME  OF  WEEK  (sec) 


FIGURE  11.  PAN  SOLUTIONS  FOR  STATIC  SITE  85401,  DAY  074 
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COMPONENT  (m) 


FIGURE  14.  PPS  SOLUTIONS  FOR  STATIC  SITE  85401,  DAY  076 


FIGURE  15.  PAN  SOLUTIONS  FOR  STATIC  SITE  85401,  DAY  076 


COMPONENT  (m) 
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NAVIGATION  SOLUTIONS  FOR  85401,  DAY  077 
SIMULATED  PPS  FROM  THE  BROADCAST  EPHEMERIDES 


GPS  TIME  OF  WEEK  (sec) 


FIGURE  16.  PPS  SOLUTIONS  FOR  STATIC  SITE  85401,  DAY  077 


NAVIGATION  SOLUTIONS  FOR  85401 ,  DAY  077 

PRECISE  ABSOLUTE  NAVIGATION 


GPS  TIME  OF  WEEK  (sec) 


FIGURE  17.  PAN  SOLUTIONS  FOR  STATIC  SITE  85401,  DAY  077 
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DYNAMIC  SITE  SOLUTIONS 

A  map  of  the  roads  at  NSWCDD  on  which  the  dynamic  data  were  collected  is  shown  in 
Figure  18.  The  two  static  site  locations  are  indicated  as  MERE  and  ASTW.  The  truth  data  for  the 
PAN  solutions  are  the  OTF  results  between  MERE  and  one  of  the  two  vehicles  that  were  part  of  the 
original  test.  The  path  the  vehicle  took  began  at  the  location  marked  Turntable,  then  followed 
McVay  Road  to  Dahlgren  Road.  At  the  Dahlgren  Road  intersection,  the  vehicles  made  a  left  turn  and 
followed  Dahlgren  Road  east  to  the  intersection  with  Caskey  Road.  A  right  turn  on  Caskey  Road 
was  followed  by  two  left  turns  on  two  secondary  roads  to  return  to  Dahlgren  Road.  After  the  Caskey 
Road  intersection,  the  vehicle  proceeded  back  to  the  starting  point  at  the  Turntable.  Five  laps  around 
this  circuit  are  shown  in  the  Figure  19  east  vs.  north  truth  plot.  A  second  truth  plot  showing  the 
ellipsoid  height  vs.  east  is  included  as  Figure  20.  This  figure  shows  that  there  is  little  change  in  height 
around  the  circuit. 

The  PPS  solutions  for  the  circuit  are  shown  in  Figures  21  and  22.  Figure  21  shows  the 
solutions  in  the  east  and  north  components,  while  Figure  22  shows  the  east  and  ellipsoid  height 
components.  The  corresponding  PAN  solutions  are  shown  in  Figures  23  and  24.  For  comparison, 
the  SPS  solutions  are  shown  in  Figures  25  and  26.  Note  the  change  in  scale  in  Figure  26.  On  the  east 
and  north  scale,  the  difference  between  the  PPS  and  PAN  solutions  are  difficult  to  appreciate,  while 
the  improvement  in  ellipsoid  height  is  easy  to  see.  The  area  marked  DETAIL  in  Figure  19  is  used  to 
illustrate  the  east  and  north  differences  on  a  smaller  scale.  For  the  detail  area,  the  east  and  north 
components  of  the  PPS  solutions  are  shown  in  Figure  27,  the  east  and  ellipsoid  heights  in  Figure  28. 
The  corresponding  PAN  solutions  are  shown  in  Figures  29  and  30.  A  plot  of  the  difference  between 
the  PPS  solutions  and  the  truth  is  shown  in  Figure  31.  Each  differenced  component  is  plotted 
separately  vs.  GPS  time  of  week.  The  corresponding  PAN  solutions  are  shown  in  Figure  32.  A  table 
summarizing  the  differences  is  included  as  Table  5. 


DIFFERENTIAL  CORRECTIONS 


Early  in  the  development  of  PAN,  it  was  thought  that  the  SPS  dynamic  vehicle  solutions 
obtained  during  the  Dahlgren  vehicle  test,  described  earlier,  could  be  made  to  simulate  the  PPS 
solutions.  These  data  were  obtained  by  four  Trimble  4000SE  receivers.  These  are  unkeyed  receivers, 
but  are  capable  of  recording  two-frequency  pseudorange  and  phase  data  even  when  AS  is  operating. 
Since  these  receivers  are  primarily  used  for  kinematic  positioning,  the  SA  effects  are  of  little 
importance  to  the  user.  In  the  case  of  navigation,  the  SA  dither  and  radial  epsilon  errors  could  be 
substantially  reduced  by  computing  range  corrections.  This  technique  was  expected  to  perform  an 
adequate  dither  correction  without  resorting  to  formal  SA  removal  in  the  DCF.  The  range 
corrections  were  to  be  computed  from  static  sites  whose  position  was  known.  These  would  then  be 
transferred  to  the  dynamic  sites  during  postprocessing. 


16 


NSWCDD/TR-95/196 


FIGURE  18.  MAP  SHOWING  ROUTE  OF  VEHICLE 
AROUND  DAHLGREN  CIRCUIT 


FIGURE  19.  TRUTH  TRACK  AROUND  VEHICLE  CIRCUIT 
CORRESPONDING  TO  FIGURE  18 
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FIGURE  20.  TRUTH  TRACK  CORRESPONDING  TO  FIGURE  18  SHOWING 
ELLIPSOID  HEIGHT  VARIATION 


FIGURE  21.  PPS  TRACK  AROUND  VEHICLE  CIRCUIT 
CORRESPONDING  TO  FIGURE  19 
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FIGURE  22.  PPS  TRACK  CORRESPONDING  TO  FIGURE  20  SHOWING 
ELLIPSOID  HEIGHT  VARIATION 


FIGURE  23.  PAN  TRACK  AROUND  VEHICLE  CIRCUIT 
CORRESPONDING  TO  FIGURE  19 
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FIGURE  24.  PAN  TRACK  CORRESPONDING  TO  FIGURE  20  SHOWING 
ELLPSOID  HEIGHT  VARIATION 


SPS  NAVIGATION  SOLUTION 

TRACK  ID  95 


FIGURE  25.  SPS  TRACK  AROUND  VEHICLE  CIRCUIT 
CORRESPONDING  TO  FIGURE  19 
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SPS  NAVIGATION  SOLUTIONS 

TRACK  ID  95 


-1100  -1000  -900  -800  -700  -600  -500  -400 

EAST  (m) 


FIGURE  26.  SPS  TRACK  CORRESPONDING  TO  FIGURE  1 9  SHOWING 
ELLIPSOID  HEIGHT  VARIATION 


FIGURE  27.  PPS  TRACK  AROUND  VEHICLE  CIRCUIT 
CORRESPONDING  TO  THE  DETAIL  IN  FIGURE  19 


21 


NSWCDD/TR-95/196 


FIGURE  28.  PPS  TRACK  CORRESPONDING  TO  DETAIL  IN  FIGURE  19 
SHOWING  ELLPSOID  HEIGHT  VARIATION 


ELLIPSOID  HEIGHT  (m) 
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PRECISE  ABSOLUTE  NAVIGATION 
TRACK  ID  95 


EAST (m) 


FIGURE  30.  PAN  TRACK  CORRESPONDING  TO  DETAIL  IN  FIGURE  19 
SHOWING  ELLIPSOID  HEIGHT  VARIATION 


PPS  SOLUTION  MINUS  TRUTH 


GPS  TIME  OF  WEEK  (S) 


FIGURE  3 1 .  DIFFERENCE  BETWEEN  PPS  SOLUTIONS  AND  TRUTH  FOR 
VEHICLE  CIRCUIT  AT  DAHLGREN 
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PAN  SOLUTION  MINUS  TRUTH 


FIGURE  32.  DIFFERENCE  BETWEEN  PAN  SOLUTION  AND 
TRUTH  FOR  VEHICLE  CIRCUIT  AT  DAHLGREN 


TABLE  5.  DffFERENCE  TABLE  FOR  VEHICLE  SOLUTIONS  AROUND  DAHLGREN  CIRCUIT 


UNIT: 

METERS 

PPS  SOLUTIONS 

PAN  SOLUTIONS 

EAST 

NORTH 

HEIGHT 

RANGE 

EAST 

NORTH 

HEIGHT 

RANGE 

MEAN 

1.73 

8.49 

-5.67 

11.20 

-0.90 

3.36 

-0.82 

4.18 

SD 

1.21 

2.59 

3.54 

1.52 

0.43 

1.17 

2.17 

1.27 

METHOD 

The  satellite  range  corrections  can  be  found  from  data  collected  at  static  sites  whose  position 
is  already  well  known.  There  were  two  such  sites  during  the  dynamic  vehicle  tests,  the  site  at 
Dahlgren  (MERE),  and  the  site  at  Corbin,  Virginia,  33  km  southwest  (ASTW).  Since  the  navigation 
solutions  are  independent  from  one  time  to  the  next,  the  range  corrections  need  to  be  obtained  in  a 
similar  manner.  If  there  are  n  satellites  in  common  view  at  site  k,  then  there  are  n  pseudorange 
observations.  The  unknowns  include  the  local  clock  bias  and  the  n  satellite  range  biases.  It  will  be 
assumed  that  the  site  coordinates  are  known.  If  the  unknown  local  clock  bias  is  lumped  in  with  the 
satellite  range  biases,  then  it  can  be  removed  from  the  problem  by  differencing. 
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One  satellite  was  selected  as  the  reference  and  arbitrarily  assigned  zero  range  bias.  The 
computed  ranges  were  constructed  from  the  satellite  ephemerides,  clock  predictions,  and  the  known 
site  position.  The  reference  satellite  range  was  subtracted  from  all  the  others  leaving  n-I 
observations  and  n-1  unknowns  that,  in  principle,  could  be  solved.  Since  there  were  actually  two 
static  sites,  there  were  2(n-l)  observations  and  n-1  unknowns.  The  problem  was  over-determined 
and  least  squares  estimation  could  be  used. 

These  range  biases  were  then  treated  as  range  corrections  and  subtracted  from  the  computed 
ranges  at  each  dynamic  site.  The  feet  that  one  satellite  was  given  a  zero  bias,  when  in  fact  it  was  not 
zero,  does  not  matter  when  the  corrections  are  made,  as  long  as  the  determination  of  the  true  vehicle 
clock  offset  is  not  important.  The  range  corrections  just  add  the  reference  satellite  range  bias  to  the 
local  clock  ofifeet,  and  the  combined  offset  is  solved  in  the  usual  way  during  the  navigation  solution. 
These  postprocessed  navigation  solutions  are  differential  solutions  because  they  depend  upon 
information  supplied  by  static  sites  whose  positions  are  known.  All  common  range  biases,  such  as 
dither,  radial  orbit  error,  residual  satellite  clock  estimation  error,  and  correlated  refractions  errors, 
are  therefore  removed  and  the  navigation  solutions  are  improved. 

RESULTS 

Examples  of  the  range  corrections  are  shown  in  Figure  33  with  PRN14  as  the  reference 
satellite.  The  degree  to  which  the  curves  are  correlated  is  the  result  of  the  common  time-varying 
range  bias  from  the  reference  satellite  and  the  local  clock. 

Figures  34  and  35  show  the  results  obtained  from  the  differentially  corrected  observations. 
These  should  be  compared  to  Figures  25  and  26,  which  show  the  original  uncorrected  SPS  solutions. 
Since  the  purpose  was  to  eliminate  only  the  dither  component  of  SA,  these  differential  solutions  could 
not  be  us^  as  the  equivalent  to  the  PPS  solutions.  The  PPS  solutions  include  the  orbit  error  from 
the  broadcast  ephemerides,  whereas  the  differentially  corrected  solutions  have  some  of  the  orbit  error 
(and  other  effects)  removed.  Table  6  summarizes  the  differences  between  the  differentially  corrected 
solutions  and  truth.  Notice  that  the  mean  errors  are  smaller  than  either  the  PPS  or  the  PAN  solutions 
listed  in  Table  5.  The  standard  deviations  are  similar  to  the  PAN  solutions. 

The  differences  between  the  SPS  and  the  differentially  corrected  SPS  solutions  vs.  time  are 
illustrated  in  Figures  36  and  37.  Notice  that  the  scale  for  the  component  axis  is  ±100  m  in  the  case 
of  the  SPS  results,  and  ±15  m  in  the  differentially  corrected  case.  Since  the  differentially  corrected 
solutions  were  not  representative  of  the  PPS  results,  they  were  not  used  for  evaluation  of  the  PAN 
solutions.  Instead,  the  data  files  and  broadcast  ephemerides  for  all  four  sites  were  corrected  by  DMA 
in  their  DCF. 
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FIGURE  33.  SATELLITE  RANGE  CORRECTIONS  DURING  DYNAMIC 
VEfflCLE  CIRCUIT  AT  DAHLGREN 


FIGURE  34.  DIFFERENTIALLY  CORRECTED  SPS  SOLUTIONS 
USING  RANGE  CORRECTIONS 
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TABLE  6.  DIFFERENCES  BETWEEN  SPS  OR  DIFFERENTIALLY  CORRECTED  SPS 

SOLUTIONS  AND  TRUTH 


UNITS: 

METERS 

SPS  SOLUTIONS 

CORRECTED  SPS  SOLUTIONS 

EAST 

NORTH 

HEIGHT 

RANGE 

EAST 

NORTH 

HEIGHT 

RANGE 

MEAN 

-5.20 

18.72 

-20.85 

56.84 

0.07 

-0.07 

0.35 

1.82 

SD 

19.65 

21.79 

50.00 

30.68 

0.62 

1.10 

1.64 

1.06 

DISCUSSION 


This  investigation  has  demonstrated  that  PPS  navigation  solutions  from  field  sites  can  be 
upgraded  to  precise-ephemerides  quality  without  reprocessing  the  original  observations.  If  files  are 
kept  of  the  unclassified  navigation  solutions,  the  time,  and  the  satelhtes  being  tracked,  and  a  source 
of  the  SA-corrected  broadcast  ephemerides  and  precise  ephemerides  are  available,  then  the  PPS 
solutions  obtained  with  the  broadcast  ephemerides  can  be  improved  at  a  later  date.  This  PAN 
capability  can  make  important  contributions  to  military  mapping  operations  where  vehicles  equipped 
with  PPS  receivers  can  gain  access  to  areas  lacking  geodetic  control  points. 
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SPS  SOLUTION  MINUS  TRUTH 


FIGURE  36,  DIFFERENCE  BETWEEN  SPS  SOLUTIONS  AND  TRUTH 
FOR  DYNAMIC  VEHICLE  95 


FIGURE  37.  DIFFERENCE  BETWEEN  CORRECTED  SPS  SOLUTIONS 
AND  TRUTH  FOR  DYNAMIC  VEHICLE  95 
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PAN  PROGRAM  DESCRIPTION 


In  developing  the  code  for  PAN,  programming  methodologies  of  portability,  readability,  and 
error  analysis  were  taken  into  consideration.  ANSI  standards  and  machine-independent  code  were 
used  to  allow  the  program  to  be  executed  on  more  than  one  platform.  Prologues  at  the  beginning  of 
each  subroutine,  comments  throughout  the  subroutines,  as  well  as  descriptive  names  for  variables  and 
routine  names  were  some  measures  taken  to  make  the  code  more  readable.  The  program  was 
compiled  and  tested  in  small  components,  debug  print  statements  were  placed  throughout  the 
program,  and  intermediate  results  were  compared  for  accuracy  to  enhance  the  error  analysis  process. 
The  code  was  also  processed  through  a  FORTRAN  code  checker  that  had  the  capabilities  to  find 
programming  errors  the  compiler  might  have  missed,  such  as  undefined  variables  and  call-line 
inaccuracies. 

The  top-down  design  approach  of  starting  with  a  problem  and  breaking  the  problem  into 
smaller  tasks  was  used  in  developing  this  code.  First,  an  outline  was  developed  fi'om  the  formulation. 
This  outline  then  was  used  as  a  basis  for  developing  pseudo-code.  Modular  routines,  where  each 
different  task  is  developed  into  its  own  subroutine,  were  developed  using  the  pseudo-code  as  a  basis. 
Routines  for  some  of  the  tasks  involving  matrix  computation  and  interpolation  already  existed  in 
mathematical  utility  libraries.  These  existing  routines  were  merged  into  a  library  directory  that  would 
later  be  linked  with  the  program.  Include  files  were  also  used  for  all  global  variables. 

PAN  has  the  capabilities  to  solve  for  one  of  the  following  options  over  a  user-specified  time 
interval: 

•  Solve  for  a  PPS  solution  using  broadcast  ephemerides, 

•  Solve  for  a  PE  solution  using  precise  ephemerides, 

•  Solve  for  both  PPS  and  PE  solutions. 

The  input  requirements  for  PAN  are: 

•  A  configuration  file 

•  1  RENEX  Navigation  file 

•  1  RINEX  Observation  file 

If  the  user  is  solving  for  PE  solution,  then  the  following  files  are  also  needed: 

•  1  to  n  Satellite  Trajectory  files  in  PC  ASCII  format 

•  1  Satellite  Clock  file  in  PC  ASCII  format 

If  the  user  is  solving  for  just  PE,  then  the  user  needs  to  input: 


A  PPS  solution  file 
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The  configuration  file  is  a  file  containing  all  user  input  parameters  needed  for  execution.  This 
file  uses  keyword  specification  in  that  it  searches  for  keywords  in  the  file  to  indicate  what  the 
parameter  values  are.  This  allows  the  user  to  put  the  parameters  in  any  order  in  the  file,  as  well  as 
to  fi'eely  put  comments  throughout  the  configuration  file.  The  keyword  search  method  also  allows 
the  user  to  omit  unnecessary  parameters  from  the  configuration  file.  With  the  ke5rword  specification 
method,  the  program  can  be  easily  modified  to  add  new  input  parameters.  This  method  also  gives 
the  user  more  flexibility  in  modifying  the  configuration  file.  A  full  description  of  the  different  input 
parameters  that  are  acceptable  by  PAN  can  be  found  in  the  PAN  User’s  Guide  section. 

The  user  can  execute  the  program  from  the  DOS  prompt  with  the  following  command: 
>pan.  exe  or  fullpath  if  not  in  directory  containing  exe 

After  the  user  types  in  the  above  command,  the  following  line  will  appear  on  the  screen,  and 
the  user  can  then  type  in  the  configuration  filename  as  indicated  by  italics; 

Enter  configuration  filename  (40  chars  max): 
c:  \pan\datapc\view.  cfg 

The  program  will  then  begin  processing  solutions  for  each  time  step  in  the  given  interval.  If 
the  user  wishes  to  have  debug  statements  printed  out,  this  option  is  available.  The  output  files  consist 
of  either  one  or  two  files,  depending  on  whether  the  option  to  compute  the  PPS  solution,  the  PE 
solution,  or  both,  was  chosen.  The  output  files  will  contain  the  X,  Y,  Z  solutions  (meters),  the  clock 
solution  (meters),  the  solution  in  longitude  and  latitude  (degrees),  the  east,  north,  and  vertical 
differences,  the  RSS  of  residuals,  and  the  GDOP  solution.  A  description  of  these  output  files  can  be 
found  in  the  PAN  User’s  Guide  section. 

PAN  was  written  in  FORTRAN  and  complied  and  executed  on  several  different  platforms 
including  several  different  PC  models,  and  an  RISC  6000.  PAN  contains  approximately  6600  lines 
of  code  separated  into  three  different  directories: 

•  DITHER  5  routines 

•  PAN  25  routines 

•  LIB  34  routines 

PAN  took  55  sec  to  complete  the  execution  of  a  case  that  was  set  up  to  solve  for  both  the  PE 
and  PPS  solution  for  six  satellites  over  a  20-sec  time  interval  with  a  step  size  of  1  sec.  In  a  similar 
case,  solving  for  both  solutions  using  25  satellites  over  a  120-sec  time  interval  using  a  time  step  of 
30  sec,  the  program  took  3  min  to  complete. 
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PAN  USER’S  GUroE 


This  program  has  the  ability  to  improve  PPS  Absolute  Navigation  Solutions  so  that  they 
approach  the  quality  that  is  achievable  from  the  DMA  precise  ephemerides  and  clock.  Real-time 
navigation  solutions  necessarily  need  to  use  the  broadcast  ephemerides  as  the  source  of  the  satellite's 
position  and  clock,  and  therefore  include  the  satellite  prediction  errors  in  the  solution.  Since  the 
precise  ephemerides  and  clock  is  a  smoothed  fit  that  is  not  available  in  real  time,  it  avoids  the 
prediction  errors  and  consequently  allows  significantly  better  results  to  be  realized. 

PROGRAM  STRUCTURE 

PAN  involves  two  major  tasks.  First  the  program  solves  the  PPS  equivalent  solution  by  the 
following  steps: 

1 .  Reads  in  RINEX  observation  files  and  satellite  ephemeris  files. 

2.  Implements  the  satellite  ephemeris  algorithm. 

3.  Chooses  the  correct  ephemeris  message  and  evaluates  the  satellite  position  at  the  time  of 
transmission  for  each  observation. 

4.  Finds  the  ionospherically  corrected  pseudorange  from  the  "raw"  pseudoranges. 

5 .  Forms  the  0^  vector  from  the  pseudorange. 

6.  Forms  the  C^  vector  from  the  broadcast  ephemeris  and  an  estimate  of  the  current  absolute 
position. 

7.  Computes  the  partial  derivatives. 

8.  Solves  for  a  navigation  solution  at  each  time  step. 

Second,  the  program  converts  the  PPS  solution  to  the  precise-ephemerides  equivalent  by  the 
following  method: 

1.  With  the  PPS  navigation  solution  at  each  time  step  and  the  broadcast  ephemerides,  it 
forms  the  computed  range  vector  Cf 

2.  With  the  PPS  navigation  solution  at  each  time  step  and  the  precise  ephemerides,  it  forms 
the  computed  range  vector  Cp. 

3.  Computes  the  partial  derivatives  from  the  PPS  solution  and  the  precise  ephemerides. 

4.  Solves  for  the  PAN  solution. 

PAN  is  set  up  so  the  user  can  either  solve  for  the  PPS  solution,  the  PPS  solution  and  PE 
solution,  or  just  the  PE  solution.  The  program  uses  the  set  of  satellite  PRN  numbers  on  the 
observation  file  to  solve  the  PPS  solution.  For  the  PE  solution,  the  user  needs  to  list  the  satellite 
PRN  numbers  to  solve  for  in  the  configuration  file.  There  must  be  a  trajectory  file  for  each  satellite 
specified,  and  these  satellite  PRN  numbers  must  appear  on  the  satellite  clock  file,  as  well  as  on  the 
Broadcast  Ephemeris  file.  If  the  user  selects  to  just  solve  for  the  PE  solution,  an  input  file  containing 
PPS  solution  for  the  time  ranch  must  be  available.  If  a  satellite  contains  events  in  the  satellite  clock 
file,  this  satellite  will  be  omitted  from  solving  the  PE  solution. 
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PROGRAM  EXECUTION 

PAN  will  ask  for  a  configuration  filename.  Enter  the  name  and  pathname  (up  to  40  charac¬ 
ters).  PAN  will  execute  over  an  interval  of  time  and  output  several  files  depending  on  the  user 
options. 


FILE  INPUT/OUTPUT 


Pan  Input 

A  configuration  file 
1  RINEX  Navigation  File 
1  RINEX  Observation  File 

1  -  MAXSAT  Satellite  Trajectory  Files  (PC  ASCII  formatted) 

(only  needed  if  solving  for  PE  solution) 

1  Satellite  Clock  File  (PC  ASCII  formatted)  (only  needed  for  PE) 

DITHER  output  file  containing  dither  estimations  (optional  file  not  needed  for  PAN) 
PSS  Solution  File  (only  needed  if  solving  for  PE  only) 


Pan  Output 

One  to  four  output  files  can  be  generated  giving  solutions  in  different  formats.  Either  a  file 
with  PPS  solution  results  or  a  file  with  PE  solution  results. 


FILE  FORMATS 


Pan  Configuration  File 

There  is  no  specific  order  to  the  placement  of  input  values  in  the  file.  The  program  searches 
for  keywords  within  the  configuration  file.  These  keywords  must  be  placed  at  the  start  of  a  line. 
Comment  lines  are  represented  by  "#"  signs  and  ignored  by  the  program. 

Here  is  a  list  of  valid  ke5rwords: 

NAV 
OBSERV 
T 
C 


RINEX  Navigation  file  name  (max  40  char) 
RINEX  Observation  file  name  (max  40  char) 
Trajectory  Filenames  (max  40  char) 

Satellite  Clock  filename  (max  40  char) 
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The  configuration  file  can  have  either  one  or  both  of  the  following  output  files. 

OTPPS  create  output  file  for  PPS  solution  (filename) 

OTPE  create  output  file  for  PE  solution  (filename) 

WEEK  GPS  week# 

TIME  Start  time  End  time  (seconds  into  GPS  week) 

RUNTYPE  Indicates  what  type  of  run  to  execute 

(1=PPS  only;  2=PPS  and  PE;  3=PE  only) 

LL  has  four  variables  on  the  line:  iopt  (which  tells  program  if  the  latitude  and  longitude  position  are 
in  degrees  or  radians),  latitude  position,  longitude  position,  and  height. 

LL  IOPT  (0=degrees  1  =radian)  LAT  LONG  HEIGHT 

The  following  are  optional  (These  can  be  put  in  the  configuration  file  if  the  user  wants  to  change  the 
default  values): 


INTERVAL 

TEMP 

PRESSURE 

HUM 

D 

PELIST 

PPSFILE 

DEBUG 


time  step  (default  value  =1) 
temperature  (default  value  =  280  KELVIN) 
pressure  (default  value  =  990  millibar) 
humidity  (default  value  =  50%) 

Dither  Estimation  File  (40  chars  max) 

List  of  satellite  PRN#s  to  solve  PE  solution  for  (default  list:  1  2  3  4  5) 
Input  PPS  filename  (40  chars  max)  (must  have  if  RUNTYPE  =  3) 

1  (option  to  turn  on  debug  statements)  (default  =  0  —  do  not  print  out) 


A  sample  PAN  configuration  is  shown  Figure  38. 


PPS  Solution  Input  File  Format 

Each  line  contains:  time  xsoln  ysoln  zsoln  clockoffset. 
format:  Ix,i7,4dl7.10 


Pan  Output  File  (for  both  PPS  and  PE  solution  files! 

One  record  for  each  time  processed  with  the  following  format: 


time  of  week  (sec) 
X  solution  (m) 
y  solution  (m) 
z  solution  (m) 
clock  solution  (sec) 
longitude  (degrees) 


— 18.0 

—  fl3.3 

—  fl3.3 

—  fl3.3 

—  fl3.10 

—  fl3.8(ex.  282.12345678) 


33 


NSWCDD/TR-95/196 


latitude  (degrees) 

—  fl3.8 

ellipsoid  height  (m) 
longitude 

—  fl2.4 

(integer  part  only) 
latitude 

—  i5  (ex.  282) 

(integer  part  only) 
longitude 

— -  i5 

(decimal  part  only) 
latitude 

—  fl2.8  (ex.  .12345678) 

(decimal  part  only) 

— fl2.8 

east  difference  (m) 

—  fl2.4 

north  difference  (m) 

—  fl2.4 

vertical  difference  (m) 

—  fl2.4 

RSS  of  residuals  (m) 

—  fl2.4 

GDOP 

—  flO.2 

#This  is  a  sample  configuration  file  using  default  values 

#  solve  for  both  PPS  and  PE  solutions 
RUNTYPE  2 

#  list  of  satellites  to  solve  PE  solution  for 
PELIST  1  14  22  25  27  29  31 

NAV  c:\pan\data\4872048b.nav 
OBSERV  c:\pan\data\4872048b.rnx 
T  c:\traj\P001736 
T  c:\traj\P002736 
T  c:\traj\P014736 
T  c:\traj\P022736 
T  c:\traj\P025736 
T  c:\traj\P029736 
T  c:\traj\P031736 
C  c:\traj\GC736 

#  This  case  is  to  include  dither  file 
D  c:\dither\dithot 

#  Want  both  output  files 
OTPPS  outl.pps 

OTPE  outl.pe 
WEEK  736 

#  time  interval  (start  and  end  time  (secs) 

TIME  411062  411065 

#  lat  long  in  radian 

LL  1  0.669d0  4.938604d0  -29.00009d0 

#End  of  sample  configuration  file 

##################################################### 


FIGURE  38.  PAN  CONFIGURATION  FILE 
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PRODUCT  DISCLAIMER 


The  mention  of  a  commercial  company  or  product  does  not  constitute  an  endorsement  by 
NSWCDD  or  DMA.  The  information  in  this  paper  is  not  authorized  for  use  for  publicity  or 
advertisements. 
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VERSION:  1.4 


PRECISE  ABSOLUTE  NAVIGATION 
FORMULATION 

Bruce  R.  Hermann  K12 

Naval  Surface  Warfare  Center,  Dahlgren  Division 
Updated:  22  September  1995 


DESCRIPTION 

This  task  will  develop  a  capability  to  improve  Precise  Positioning  Service  (PPS)  absolute  navigation 
solutions  so  that  they  approach  the  quality  that  is  achievable  from  the  Defense  Mapping  Agency 
(DMA)  predse  ephemerides  and  clock.  Real-time  navigation  solutions  necessarily  use  the  broadcast 
ephemerides  as  the  source  of  the  satellite's  position  and  clock  and,  therefore  include  satellite 
prediction  errors  in  the  solution.  Since  the  precise  ephemerides-and-clock  is  a  smoothed  fit  that  is 
not  available  in  real  time,  it  avoids  prediction  errors  and  consequently  allows  significantly  better 
results  to  be  realized. 

ASSUMPTIONS 

The  Qobal  Positioning  System  (GPS)  time  and  the  PPS  solutions  to  be  improved  are  given,  as  is  the 
receiver  type  and  the  version  of  the  navigation  algorithm  used.  Algorithm  information  will  be  used 
to  establish  which  satellites,  of  those  in  view,  produced  the  particular  navigation  solution.  To 
compute  the  difference  between  the  ephemerides,  an  independent  source  of  the  Selective  Availability- 
(SA-)corrected  broadcast  ephemerides  and  the  DMA  precise-ephemerides-and-clock  are  required. 
Finally  for  testing  purposes,  raw  pseudorange  and  phase  observations  from  a  typical  receiver  will  be 
needed  in  addition  to  the  information  above. 


METHOD 

Satellite  ephemerides,  clocks,  and  the  PPS  solution  are  given.  The  form  of  the  observation  equations 
that  produced  the  navigation  solution  can  be  reconstmcted  as  in  Equation  (A-1).  The  original  vector 
of  observation  was  Og,  the  computed  ranges  were  in  the  Cg  vector,  and  the  initial  state  vector  was 
Xg.  None  of  these  vectors  were  saved. 

SC 

=  Co  ^  ^  LXo  (A-l) 

dXg 


The  solution  toEquation  (A-l)  is  AXg,  from  which  the  new  state  vector  Aj  (the  PPS  solution)  is 
foimd:  Xf^Xg+AXg.  Since  the  same  observations  would  be  used  regardless  of  the  ephemerides,  the 
observation  equations  using  the  broadcast  ephemerides  (subscript  b)  can  be  equated  to  the  equations 
using  the  precise  ephemerides  (subscript  p)  inEquation  (A-2). 


dX. 


AX.  =  C  +  — i 
‘  ^  dX. 


AX 


(A-2) 


The  given  PPS  navigation  solution  Xg  is  used  in  the  computation  for  the  vector  C  on  both  sides  of 
Equation  (A-2).  Assuming  operation  in  the  linear  region,  the  solution  from  Equation  (A-l)  used 
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VERSION:  1.4 


along  with  the  original  satellite  ephemerides  to  compute  Q  should  require  no  further  correction.  This 
makes  AAj  »  0.  Therefore  Equation  (A-2)  can  be  rearranged  into  a  form  similar  to  that  in 
Equation  (A-1).  The  result,  Equation  (A^3),  can  be  solved  in  the  usual  least  squares  sense  to  obtain 
the  correction  AXp  that  should  be  applied  to  the  given  to  form  the  precise  solution:  Xp=Xi+AXp. 


C.  -  C=  — 


p 


dXp  ^ 


(A-3) 


TASK  OUTLINE 

PPS  Equivalent  Solution 

1 .  Read  RINEX  observation  files  and  satellite  ephemeris  files. 

2.  Implement  the  satellite  ephemeris  algorithm. 

3.  Choose  the  correct  ephemeris  message,  and  evaluate  the  satellite  position  at  the  time  of 
transmission  for  each  observation. 

4.  Check  the  phase  observations  for  cycle  jumps. 

5.  Smooth  the  pseudoranges  with  the  phase  observations. 

6.  Form  the  Oq  vector  from  the  smoothed  pseudoranges  for  use  in  Equation  (A-1). 

7.  Form  the  Cg  vector  from  the  broadcast  ephemeris  and  an  estimate  of  the  current  absolute 
position. 

8.  Compute  the  partial  derivatives  for  use  in  Equation  (A-1). 

9.  Solve  for  a  navigation  solution  at  each  time  step  and  save  the  results. 

Convert  PPS  Solution  to  the  Precise  Ephemerides  Equivalent: 

1 .  With  the  PPS  navigation  solution  at  each  time  step  and  the  broadcast  ephemerides,  form  the 
computed  range  vector,  Q  in  Equation  (A-3). 

2.  With  the  PPS  navigation  solution  at  each  time  step  and  the  precise  ephemerides,  form  the 
computed  range  vector,  Cp  in  Equation  (A-3). 

3.  Compute  the  partial  derivatives  from  the  PPS  solution  and  the  precise  ephemerides. 

4.  Solve  Equation  (A-3)  for  the  precise  absolute  navigation  solution. 

5.  Compare  the  PPS  solutions  with  the  truth. 

The  code  should  be  developed  in  FORTRAN  or  C  so  that  it  can  be  compiled  on  a  single-user  PC  as 
well  as  a  workstation.  The  data  files  for  the  equivalent  PPS  solution  are  currently  in  RINEX  format 
on  PC-compatible  disks.  The  precise  ephemerides  are  available  from  OMNIS  files  for  GPS  week 
736. 

CONFIGURATION  INPUT  FILE 

The  user  will  interface  with  the  program  through  a  configuration  file.  This  ASCII  file  will  contain 
the  information  needed  by  the  program  to  perform  the  computations.  The  user  can  keep  copies  of 
several  files  with  different  names  for  different  data  spans.  The  file  will  include  the  following 
information: 

1 .  RINEX  data  file  name 

2.  RINEX  navigation  file  name. 

3.  Precise  ephemeris  file  names,  one  per  satellite. 
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4.  Precise  clock  file  name. 

5.  Initial  site  position;  longitude,  latitude,  height  (deg,  deg,  m). 

6.  Start  and  stop  times  (GPS  time  of  week  in  seconds). 

7.  A  list  of  satellite  PRN  numbers  to  omit  from  processing. 

GENERATE  AN  EQUIVALENT  PPS  NAVIGATION  SOLUTION 

In  the  operational  environment,  the  PPS  navigation  solutions  will  be  obtained  from  the  field. 
Presumably  these  solutions  were  obtained  with  a  receiver  that  can  correct  for  the  effects  of  Selective 
Availability  (SA)  and  AntiSpoofing  (AS).  AS  encrypts  and  SA  corrupts  the  GPS  signals  making  the 
precise  signal  and  the  full  accuracy  of  the  broadcast  ephemerides  accessible  to  authorized  users  only. 
These  receivers  can  display  the  navigation  solution  in  real  time,  but  the  field-corrected  satellite 
ephemerides  and  observations  are  classified  and  therefore  are  not  available  to  the  user.  For  this 
demonstration,  the  equivalent  to  a  PPS  solution  needs  to  be  computed  fi-om  known  observations,  with 
a  known  algorithm,  and  with  known  satellite  ephemerides. 

The  equivalent  PPS  solution  used  to  demonstrate  the  method  is  derived  from  data  that  do  not  have 
the  effects  of  SA  (epsilon)  removed  fi-om  the  broadcast  ephemerides  or  SA  (dither)  from  the  satellite 
fi-equency  standard,  but  does  provide  dual  frequency  pseudorange.  These  data  were  obtained  with 
Trimble  4000SSE  receivers  as  part  of  the  Evans  Drive-By  experiment  performed  in  February  1994."^'^ 
The  Trimble  4000SSE  receiver  gives  Lj  C/A  and  the  cross  correlation  between  Lj  and  Lj  Y  codes 
(when  AS  is  implemented)  to  provide  the  L2  pseudorange  observation.  The  dual  frequency  phase 
observations  are  also  available. 


GPS  OBSERVATION  MODEL 

For  best  results,  on-the-fly  (OTF)  relative  positioning  requires  the  use  of  both  pseudorange  and  phase 
observations  at  Lj  and  £-2  frequencies.  An  observation  model  that  describes  the  pseudorange  was 
presented  by  Braasch.^'^  See  Figure  A-1 .  The  notation  will  be  modified  for  the  purposes  here,  but 
the  general  form  will  be  retained.  For  each  parameter,  a  dependance  upon  satellite  will  be  expressed 
as  a  superscript  j  or  k,  a  dependance  upon  receiver  will  be  expressed  by  a  subscript  m  or  n.  The 
transmission  frequency  is  denoted  by  the  subscript /,  which  may  be  either  a  1  or  2,  and  the  time  of 
the  observation  by  the  subscript  i.  Absence  of  a  dependance  is  represented  by  a  dot. 


Braasch's  observation  equation  includes  the  following  terms;  the  true  geometric  range  r’„  i  between 
receiver  antenna  m  and  the  satellite  j  at  epoch  /, ;  the  receiver  time  offset  , ;  the  satellite  time  offset 
t'  ,  ;  propagation  effects  due  to  the  troposphere  ,  and  the  ionosphere  =  (X^f/c)  ;  the  User 

Range  Error  URE;  multipath  error  receiver  satellite  channel  delay  receiver  measurement  bias 
errors  4^4*  receiver  noise;  and  the  deliberately  introduced  effects  of  Selective  Availability  SA/  j .  In 
the  following  development,  the  User  Range  Error,  multipath  error,  receiver  channel  delay,  bias  errors, 
and  noise  wifl  be  omitted,  but  a  term  is  added  to  account  for  relativistic  effects  The  ionospheric 
constant  ,  is  related  to  the  columnar  electron  content.  The  pseudorange  observation  for 
satellite  j,  receiver  m,  frequency  1,  and  observation  i,  with  meters  as  units,  is  written  in 
Equation  (A-4). 


pL/  =  vi.,  * 


<,)  *  Ti,  *  ^  SA 


(A-4) 


A-5 


NSWCDD/TR-95/196 


VERSION:  1.4 


FIGURE  A-1 .  ONE-SATELLITE,  TWO-SITE  GEOMETRY 


A  similar  expression  can  be  used  for  the  phase  observation  if  an  integer  cycle  ambiguity  term  N„,o 
is  added.  It  represents  the  number  of  full  cycle  counts  recorded  by  the  receiver  on  a  particular 
channel  and  is  arbitrary.  Since  the  integer  cycle  term  is  part  of  the  observation,  it  will  be  included  on 
the  left-hand  side  in  Equation  (A-5).  The  phase  (cycles)  represents  the  observed  range 
(meters)  divided  by  the  wavelength  (meters),  so  the  result  has  units  of  cycles.  BothiV'„;o  and 
'have  units  of  cycles,  so  the  wavelength  A,/  is  needed  to  change  these  cycles  back  to  meters.  Note  that 
the  ionospheric  refraction  in  Equation  (A-5)  has  the  opposite  sign  than  it  had  in  Equation  (A-4).'''^ 


*  flLl 


tL  -  —kL 


+  SA 


(A-5) 


The  geometric  range  ,  is  the  magnitude  of  the  vector  to  the  satellite,  in  the  WGS84  coordinate 
system,  minus  the  vector  to  the  receiver  antenna,  as  shown  in  Equation  (A-6).  The  satellite 
coordinates  as  a  function  of  time  are  available  from  the  satellite  ephemerides,  while  the  receiver 
antenna  coordinates  are  to  be  found. 


'ii  -  V.  =  ^  -  yn,.i)y  + 


(A-6) 


SMOOTHING  PSEUDORANGE  WITH  PHASE  (Not  Implemented) 

As  shown  in  Equation  (A-5),  the  phase  observations  have  an  unknown  tracker  bias  N'„,o  that  is  site-, 
satellite-,  and  frequency-dependent.  Therefore,  the  phase  observations  cannot  be  used  to  obtain  a 
navigation  solution  in  the  same  manner  as  the  pseudoranges,  because  there  is  an  added  unknown 
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range  bias  introduced  by  each  satellite.  This  is  unfortunate,  because  the  noise  associated  with  a 
single-frequency  pseudorange  observation  is  about  an  order-of-magnitude  greater  than  the  noise  on 
the  equivalent  phase  observation.  Multipath  is  also  less  pronounced  with  phase.  For  these  reasons, 
when  pseudoranges  are  needed  but  an  instantaneous  real-time  navigation  solution  is  not  the  primary 
objective,  it  is  advantageous  to  use  the  phases  to  smooth  the  pseudoranges. 

Single  Frequency  Observations 

In  order  to  smooth  the  pseudoranges,  continuous  uninterrupted  phase  observations  from  all  satellites 
are  required  over  the  entire  period  of  interest.  The  idea  is  to  reduce  the  error  on  the  initial 
pseudorange  observation  by  averaging.  In  effect,  the  pseudoranges  at  any  time  become  the  sum  of 
the  initial  averaged  pseudorange,  plus  the  difference  between  the  phase  observation  at  and  the  initial 
phase  at  tg.  In  practice,  phase  observations  will  be  interrupted  from  time  to  time  for  a  variety  of 
reasons;  consequently,  it  is  prudent  to  plan  for  these  interruptions.  To  do  this,  data  from  each 
satellite  must  be  processed  independently,  because  the  phase  interruptions  can  occur  at  different  times 
on  different  satellites.  Over  the  time  period  of  interest,  several  interruptions  can  be  expected  to  occur 
on  a  given  satellite;  so  the  pseudorange  averaging  must  be  recomputed  after  each  continuous  span. 
The  smoothing  process  is  demonstrated  in  the  equations  that  follow. 


By  differencing  between  Equations  (A-4)  and  (A-5),  all  the  terms  on  the  right-hand  side  are  removed 
except  for  the  ionospheric  component.  This  result  for  i,  is  shown  in  Equation  (A-7). 


(A-7) 


By  differencing  in  time  two  equations  like  Equation  (A-7),  the  unknown  integer  N„,o  can  be  removed. 
Equation  (A-8)  shows  the  difference  between  observations  at  i,  and  tg. 


pL  -  -  xi,) 


(A-8) 


Equation  (A-8)  can  be  rearranged  so  that  all  the  terms  relating  to  tg  are  on  the  left-hand  side.  They 
are  expressed  in  terms  of  the  pseudorange  and  ionospheric  term  at  each  succeeding  time,  plus  the 
phase  differences  as  shown  in  Equation  (A-9). 


The  average  pseudorange,  denoted  by  the  overbar  in  Equation  (A- 10),  is  evaluated  at  tg  by  summing 
Equation  (A-9)  over  all  K  observation  times. 
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Now  the  pseudorange  at  each  time,  ^  can  be  reconstructed  by  rearranging  Equation  (A-9)  and 
inserting  the  average  value  for  the  pseudorange  at  tg  from  Equation  (A-10).  The  result, 
Equation  (A-1 1),  is  the  smoothed  pseudorange  observation  equation  for  Lj.  If  Lj  observations  are 
available,  the  same  procedure  can  be  applied  to  them.  In  cases  where  two-frequency  data  are 
available,  it  may  be  useful  to  remove  the  ionospheric  term.  This  procedure  is  discussed  in  the  next 
section. 


(A-11) 


Dual  Frequency  Observations 

When  dual-frequency  observations  are  available,  the  ionospheric  term  may  be  removed  before 
computing  the  smoothed  pseudorange.  The  pseudorange  time  difference,  derived  from 
Equation  (A-4),  is  shown  in  Equation  (A- 12);  and  the  phase  time  difference,  derived  from 
Equation  (A-5),  is  shown  in  Equation  (A-13).  Assuming  that  there  are  no  phase  jumps  between 
epochs  tg  and  r, ,  the  integer  cycle  term  N'^j.  subtracts  out  in  Equation  (A-13). 


pL,  -  pLo  =  -iKL  -Kig)^ 

c 

Lo)  +  (vi,  -  vio)  +  -  T^,)  -  c(x^.,  - 

(.tL  -  tLo)  *  (.SAi  -  SA^g) 


(A-12) 
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The  ionospheric  term  can  be  computed  from  Equation  (A-13)  if  two-frequency  phase  data  are 
available.  First,  form  the  differences  as  in  Equation  (A-13)  for  both  Lj  and  L2.  Then  difference  the 
two  time  differences.  All  the  terms  on  the  right  in  Equation  (A-13)  are  independent  of  frequency 
except  the  ionospheric  term.  The  result  of  the  difference  is  written  as  Equation  (A-14). 
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Each  of  the  ionospheric  terms  is  modeled  by  a  frequency-independent  parameter  which 
represents  the  electron  content  along  the  line  of  sight,  between  each  receiver  antenna  and  each 
satellite.  A  solution  for  the  constants  is  possible  in  terms  of  the  phase  observations  and  their 
respective  wavelengths  as  shown  in  Equation  (A- 15). 
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If  Equation  (A-15)  is  multiplied  by  2(X^,  /c)  and  added  to  both  sides  of  Equation  (A- 13),  a  time- 
differenced  observation  is  created  that  has  the  same  right-hand  side  as  the  pseudorange  time 
difference  Equation  (A-12).  Equation  (A-16)  shows  that  this  form  has  the  same  functional  properties 
as  the  pseudorange,  but  with  the  low  noise  properties  of  phase.'^'^ 
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Since  the  right-hand  side  of  Equation  (A-16)  is  the  same  as  the  right-hand  side  of  Equation  (A-12), 
the  left-hand  sides  of  those  two  equations  must  be  equal.  Exploiting  this,  the  pseudorange  at  /,  can 
be  rewritten  in  terms  of  the  two-frequency  phase  observations  and  the  pseudorange  at  tg  as  in 
Equation  (A- 17).  A  similar  expression  can  be  written  for  Lj. 


pL  =  PiiO  +  *  2X5 


kl  -  k] 


(A-17) 


Equation  (A-17)  can  be  rearranged  to  find  an  average  value  for  p’^jo,  as  was  done  for  the  single 
frequency  case.  This  average  can  then  be  substituted  into  Equation  (A-17)  and  the  smoothed 
pseudoranges  computed  from  the  phases.  Equation  (A- 18)  is  used  to  compute  this  average 
pseudorange  at  to- 
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Though  it  would  be  possible  to  use  Equation  (A-4)  and  its  L2  equivalent  to  find  the  PPS  solution  for 
the  demonstration,  smoothing  with  the  phase  observations,  as  in  Equation  (A- 18),  is  a  better 
approach.  The  only  problem  is  that  phase  jumps  need  to  be  found  and  appropriate  action  taken. 

Interruptions  in  the  phase  observations  may  be  identified  by  cycle  slips  or  by  actual  time  gaps  in  the 
observations  for  a  given  satellite.  Time  gaps  in  the  observations  will  be  obvious,  but  the  identification 
of  cycle  slips  is  more  diflBcult.  In  order  to  detect  cycle  slips  in  the  phase  data  and  also  to  improve  the 
pseudorange  observations,  the  observations  between  time  epochs  can  be  differenced.  A  cycle  slip  is 
a  phenomenon  that  occurs  when  the  receiver  mistracks  the  satellite  signal  without  losing  lock  (no  time 
gap)  and  adds  or  subtracts  some  multiple  of  an  integer  (or  half  integer,  depending  upon  the  receiver 
design)  to  the  phase  observation.  Since  this  occurrence  adds  or  subtracts  a  wavelength  (or  more)  of 
range  between  the  receiver  antenna  and  the  satellite,  an  erroneous  position  solution  will  result  if  it 
remains  undetected. 

Phase  Jump  Detection 

Occasionally  the  receiver  may  slip  or  jump  cycles  when  recording  the  phase  observations  while 
tracking  a  satellite.  If  this  occurs,  the  resulting  average  pseudorange,  computed  from 
Equation  (A- 18),  will  be  in  error.  In  order  to  prevent  this  problem,  the  raw  phase  observations  need 
to  be  checked  for  consistency.  This  can  be  accomplished  through  the  use  of  Equation  (A- 14).  The 
expression  on  the  left  of  the  equals  sign  in  Equation  (A- 14)  differences  the  difference  between  the 
phase  at  the  current  time  and  the  phase  at  a  given  epoch  between  the  two  frequencies.  Each  phase 
difference  is  converted  to  units  of  meters  by  multiplying  by  the  appropriate  wavelength.  As  the  right- 
hand  side  of  the  equation  shows,  this  is  a  measure  of  the  change  in  the  ionospheric  path  length  from 
site  m  to  satellite  j  from  the  epoch  time  to  the  current  time  as  the  satellite  moves  across  the  sky. 
Generally,  this  change  is  very  smooth;  in  fact,  if  there  were  no  ionosphere,  the  right-hand  side  of  the 
equation  would  be  zero. 

A  single  cycle  phase  jump  in  or  shows  up  as  a  wavelength  step  in  the  smooth  change  due 

to  the  ionosphere.  This  jump  should  be  easily  detectable  if  it  occurs.  Since  either  or  both  phases  may 
jump  at  the  same  time,  the  minimum  jump  magnitudes  include  the  following  possibilities:  ±Ai,  iAj, 
AjiAz,  or  AftAj.  Figure  A-2  shows  an  example  of  Equation  (A-14)  for  a  short  span  of  observations 
with  no  jumps  recorded.  Figure  A-3  shows  the  same  data  span  with  several  single  wavelength  jumps 
artificially  included  in  the  observations.  PRN03  illustrates  a  positive  Aj  jump.  PRN14  illustrates  a 
positive  A2  jump.  PRN18  illustrates  a  jump  of  the  sum  Aj  +  A2;  PRN22  illustrates  a  jump  of  the 
difference  A2  -  Aj;  and  PRN25  illustrates  a  jump  of  the  difference  A,  -  A2.  PRN28  and  PRN29 
exhibit  no  change.  The  least  obvious  jump  is  the  jump  with  the  smallest  magnitude,  the  Aj  +  A2  case. 
This  is  about  5  cm,  which  is  equal  to  the  difference  in  wavelength  (24-19  cm).  Usually,  even  this 
small  jump  should  be  detectable. 

The  procedure  proposed  to  deal  with  these  jumps  is  as  follows.  Fit  a  polynomial  to  the  computation 
on  the  left-hand  side  of  Equation  (A-14)  for  each  satellite.  The  degree  of  the  polynomial  required 
will  depend  upon  the  time  interval  between  observations  and  the  dynamics  of  the  geometry  change 
between  the  satellite  and  the  site.  For  most  cases,  a  linear  polynomial  of  the  form  illustrated  in 
Equation  (A- 19)  will  suffice. 
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ACCUMULATED  IONOSPHERIC  PHASE  DIFFERENCE 
EXAMPLE  OF  NO  PHASE  JUMPS 


GPS  TIME  OF  DAY  (sec) 


FIGURE  A-2.  AN  EXAMPLE  OF  EQUATION  (A- 1 4)  WITH  NO  JUMPS 


ACCUMULATED  IONOSPHERIC  PHASE  DIFFERENCE 
EXAMPLE  OF  SEVERAL  PHASE  JUMPS 


GPS  TIME  OF  DAY  (sec) 


FIGURE  A-3.  ARTIFICIAL  PHASE  JUMPS  INSERTED  INTO  THE  DATA 
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A  tninimnm  of  two  observations  is  needed  to  evaluate  the  two  coefficients.  Perform  a  sliding  fit  over 
the  last  300  sec  of  data  or  at  least  two  observations,  whichever  is  longer.  If  a  jump  occurred  at  the 
current  time,  it  should  be  readily  evident  by  evaluating  the  residuals  of  the  fit  (the  difference 
between  the  polynomial  evaluated  at  and  the  phase  difference  observation  at  that  time).  If  this 
residual  exceeds  the  ±5  cm  minimum,  then  the  pseudorange  smoothing  must  be  restarted  for  this 
satellite.  That  means  the  sum  in  Equation  (A-18)  is  ended  and  all  the  pseudoranges  up  to  the  current 
time,  evaluated  using  Equation  (A- 17).  The  reference  time  tg  is  then  reset  to  the  current  time  and  a 
new  sum  begun. 

COMPUTE  SATELLITE  POSITION 

The  satellite  position  will  ultimately  be  computed  with  both  the  broadcast  ephemerides  and  the  DMA 
precise  ephemerides.  Both  are  Earth-Centered,  Earth-Fixed  (ECEF)  and  need  to  be  evaluated  to 
obtain  the  satellite  position  at  the  time  of  transmission.  The  procedure  for  doing  this  is  described  in 
this  section. 

Broadcast  Ephemeris 

The  broadcast  ephemerides  from  each  satellite  in  view  are  available  in  the  RINEX  file  format.  A 
satellite  position  at  the  time  of  transmission  is  required  in  order  to  compute  the  range  along  the  line 
from  the  receiving  antenna  to  the  transmitting  antenna.  Each  observed  pseudorange  and  phase 
measurement  is  tagged  with  the  local  time  of  reception  determined  by  each  receiver.  All  received 
epochs  fi-om  different  sites  at  the  same  time  will  have  the  same  time  tag,  but  the  true  time  of  reception 
at  each  site  will  be  in  error  by  an  unknown  receiver  offset.  The  times  of  transmission  from  a  given 
satellite  will  be  different  for  each  receiver  since  each  is  at  a  different  range  from  that  satellite.  The 
true  time  of  transmission  for  each  satellite  at  each  receiver  must  be  found  through  an  iterative 
process.  This  will  be  described  later. 

Once  the  time  of  transmission  is  known,  the  satellite's  ECEF  position  is  found  by  evaluating  a  series 
of  equations  with  coefficients  obtained  fi-om  the  current  ephemeris.  The  equations  are  repeated 
below.''"‘ 

WGS84  Constants 

Universal  gravitational  constant  times  the  mass  of  the  Earth:  p  =  3. 986005  xlO*"*  mVs^ 

Earth's  rotation  rate:Qe  =  7.2921151467x10"*  rad/s 
Speed  of  light  in  vacuum:  c  =  299792458  m/s 
Relativistic  constant:  -4.442809305x10"^°  s/m'^' 
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Ephemeris  Parameters  Included  in  the  Broadcast  Message 


Mean  anomaly  at  the  reference  time 
Mean  motion  difference 
Eccentricity 

Semimajor  axis  square  root 
Right  ascension  at  reference  time 
Inclination  at  reference  time 
Argument  of  perigee 

Rate  of  right  ascension 
Rate  of  inclination 

Cosine  harmonic  for  argument  of  latitude 
Sine  harmonic  for  argument  of  latitude 
Cosine  harmonic  for  orbit  radius 
Sine  harmonic  for  orbit  radius 
Cosine  harmonic  for  inclination 
Sine  harmonic  for  inclination 
Reference  time  ephemeris 


Mo 

An 


h 

6> 

Q 

^dot 

c 


a 

c„ 


Equations  for  Computation  of  the  Satellite  Position 


Semimajor  axis  of  the  Earth 
Computed  mean  motion 


a  =  (Va)^ 


'i 


n  =  +  Ln 


Mo* 


Time  from  ephemeris  reference  epoch 
Corrected  mean  motion 
Mean  anomaly 
Eccentric  anomaly 

Since  the  eccentric  anomaly  E^.  appears  on  left  of  the  equals  sign  in  the  equation  above  and  also  on 
the  ri^t  as  the  argument  of  a  sine  ftinction,  it  must  be  solved  by  iteration.  Start  with  »  M^.  on  the 
right  and  compute  a  new  E^,  Continue  iterating  until  the  new  E^  differs  from  the  previous  by  less  than 
10'^  radians  or  four  iterations,  whichever  comes  first. 


+  e  sin 


True  anomaly 

Argument  of  latitude 
Argument  of  latitude  correction 
Radius  correction 
Inclination  correction 
Corrected  argument  of  latitude 
Corrected  radius 
Corrected  inclination 


tan 


\j\  -  e^sin 


cos  jEj  -  e 


V,  +  CD 


=  C^sin  2(|>j  +  C^cos  2<j)j 

C.^cos  2(|>j  +  C^sin  2^^ 

rj=a(l-ecos  E^,)  +  br^^ 
h  =  'o  *  ^h*  ^doA 
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Positions  in  the  orbital  plane 
Corrected  longitude  of  ascending  node 

Earth  fixed  coordinates 

Satellite  clock  correction 
Relativistic  correction 


Oj  =  Op  +  (Q  - 

=  ^fccos  Q,  -  y^cos  i^sm 
=  ^fcSin  Qj  +  y,cos  /jcos 

z,  =  y^sin  i, 

=  "o  +  *  °2(*  -  *oc)^ 

Lti  =  Fe  sfa  sin 


Test  the  values  of  t-  and  /-  in  the  following  way.  If  either  difference  is  greater  than  302,400 
sec,  subtract  604,800.  If  either  difference  is  less  than  -302,400,  add  604,800. 

Precise  Ephemeris 

The  precise  satellite  ephemerides  and  clock  corrections  are  available  in  an  OMNIS  file  formaf^'®  for 
each  15 -min  epoch  during  the  8-day  fit  span.  Each  satellite  ephemeris  file  contains  the  GPS  time  of 
week  and  the  ECEF  satellite  position  and  velocity  at  that  time.  To  evaluate  the  satellite  position  at 
other  times,  an  eight-point  Lagrange  interpolation  routine  is  normally  used.  The  separate  clock  file 
contains  the  satellite  clock  time  offsets  for  all  the  satellites  at  15-niin  epochs.  These  time  offsets  can 
be  used  directly  (after  changing  the  units  to  seconds)  for  ,  in  Equation  (A-4). 

The  second  record  of  the  clock  file  contains  a  list  of  satellite  jumps  that  occur  during  the  period 
covered  by  the  ephemerides.  Since  an  interpolation  routine  will  give  erroneous  values  when  it  tries 
to  interpolate  across  a  jump,  a  satellite  with  a  jump  within  the  time  interval  to  be  processed  cannot 
be  used.  Therefore,  before  beginning  any  processing,  check  the  precise  clock  file  for  satellite  clock 
jumps  between  the  start  time  -  tj)  and  stop  times  (/„  +  tj)  listed  on  the  configuration  file.  The  value 
of  tj  is  the  extra  number  of  seconds  needed  to  accommodate  the  width  of  the  interpolating  function. 
If  such  jumps  exist,  print  a  message  to  the  user  and  then  add  that  PRN  number  to  the  list  of  satellites 
to  be  omitted  (also  listed  on  the  configuration  file). 

The  DMA  precise  ephemerides  gives  the  location  of  the  center  of  mass  of  the  GPS  satellites,  whereas 
the  observations  are,  by  definition,  referenced  to  the  phase  center  of  the  transmitting  antenna.  The 
phase  center  and  the  center  of  mass  are  different  points,  and  a  correction  must  be  added  to  account 
for  this  difference.  No  corresponding  correction  is  needed  for  the  broadcast  ephemerides  because 
it  is  already  referenced  to  the  satellite  antenna  phase  center.  An  algorithm  developed  to  compute 
this  correction  follows. 

First,  the  position  vectors  of  the  sun  and  the  satellite  need  to  be  known  at  the  time  of  interest.  The 
satellite  position  vector  is  obtained  from  the  precise  ephemeris  in  the  usual  way,  but  the  position 
vector  of  the  sun  must  be  converted  from  inertial  space.  Two  subroutines  to  perform  this  task  are 
SUNPOS  and  HAA.  The  first  computes  the  sun's  position  in  inertial  space,  and  the  second 
transforms  it  into  ECEF  coordinates.  Given  these  two  vectors:  the  ECEF  vector  for  the  sun,  and 
r  j,  the  ECEF  vector  for  the  satelhte,  the  unit  vectors  for  a  body-centered  coordinate  system  can  be 
computed.  It  is  known  that  the  GPS  antennas  always  must  point  toward  the  earth.  This  axis  is  the 
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z  body  axis  and  is  coincident  with,  but  in  the  opposite  direction  to,  the  satellite  vector  r  J.  It  is 
also  known  that  the  solar  panels,  they  body  axis,  must  be  perpendicular  to  the  earth-sun  vector.  The 
X  body  axis  completes  a  right  handed  system.  The  relationships  are  shown  in  the  following  set  of 
equations. 


r^. 


ys 


i-j  ^ 


(A-20) 


These  body-centered  unit  vectors  are  expressed  in  terms  of  the  ECEF  coordinates  in  Equation  (A-20). 
For  Block  n  satellites,  the  antenna  phase  center  lies  at  (0.279, 0,  1 .943 }  in  the  body-centered  system. 
Therefore  the  vector  to  be  added  to  the  satellite  vector  r J  (the  range  decreases)  is  given  in 
Equation  (A-21)  using  the  unit  vectors  given  in  Equation  (A-20). 

~  (•^■21) 


Satellite  Position  At  The  Time  Of  Transmission 

Let  tg  be  the  GPS  time  of  reception  (seconds  of  the  week  at  epoch  /)  as  read  from  the  RINEX  file 
containing  the  data  observed  at  a  given  site.  The  ECEF  position  of  the  satellite  at  the  time  of 
transmission  can  be  computed  using  the  equations  listed  above.  In  order  to  find  the  propagation  path 
length,  the  position  of  the  satellite  at  the  time  of  transmission  tj-  is  required.  The  expression  defining 
the  propagation  path  from  the  user  located  at  ,  to  the  satellite  at  ,  is  given  in  Equation  (A-22). 
This  is  Equation  (A-6)  with  the  GPS  times  of  signal  transmission  tj-  and  signal  reception  4  written 
explicitly. 

(A-22) 


Since  the  observation  time  tags  are  at  4,  this  frame  of  reference  will  be  chosen  as  the  earth-fixed 
frame  at  4.  This  allows  the  station  position  vector  to  be  fixed  and  defined  as  shown  in 
Equation  (A-23). 

Vn-ty  *  (A-23) 


The  time  of  transmission  is  unique  for  each  satellite  (subscript  j)  and  must  be  computed  starting  from 
the  time  of  reception.  The  GPS  time  of  transmission  is  found  by  subtracting  the  propagation  path 
delay  from  the  time  of  reception  and  including  the  correction  for  the  satellite  clock  offset.  This  is 
shown  in  Equation  (A-24). 
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f  J.  ~  fjj  - 


(A-24) 


Since  the  time  of  transmission  appears  on  both  sides  of  the  equation  (see  Equation  (A-20)),  some 
iteration  will  be  required.  For  the  first  trial,  the  propagation  term  can  be  approximated  by 

C 

0.075  sec,  and  Af/,  the  satellite  clock  correction,  by  zero.  This  will  give  a  transmission  time  to  use 
in  the  ephemeris  computations  and  allow  a  satellite  vector  and  a  clock  correction  to  be  found.  Since 
this  is  an  ECEF  computation,  the  satellite  position  is  computed  in  the  fi'ame  at  the  time  of 
transmission.  This  result  needs  to  be  rotated  into  the  fi'ame  at  the  time  of  reception  so  that  both 
vectors  are  known  in  a  common  frame.  Then  the  subtraction  shown  in  Equation  (A-22)  can  be 
performed. 


x/  =  Xycos(0)  +  y/sin(0) 
Vr  =  -  x/sm(0)  +  y^cos(Q) 


(A-25) 


The  transformation,  shown  in  Equation  (A-25),  just  rotates  the  ECEF  frame  around  the  z  axis  to 
compensate  fijr  the  rotation  of  the  earth  during  the  time  it  took  the  signal  to  travel  fi’om  the  particular 
satellite  J  to  the  receiver  antenna.  The  value  of  the  angle  is  shown  in  Equation  (A-26).  Continue  to 
iterate  Equations  (A-24)  and  (A-25),  by  computing  a  new  satellite  vector  and  clock  offset,  until  the 
change  in  0  from  one  iteration  to  the  next  is  less  than  a  microradian. 

e  =  d.  (A-26) 

c 


=  7.2921151467  xlO'*  rad  Is  (A-27) 


TROPOSPHERIC  REFRACTION  MODEL 

Tropospheric  refraction  is  caused  by  the  variation  of  the  index  of  refraction  from  that  of  fi-ee  space 
as  a  result  of  the  molecules  in  the  atmosphere.  This  causes  a  delay  and  a  bending  of  the  ray  path  as 
it  propagates  through  the  atmosphere  along  the  line  of  sight  from  the  satellite  to  the  receiving 
antenna.  This  refraction  is  independent  of  frequency  and  cannot  be  eliminated  by  the  dual  frequencies 
transmitted  by  the  GPS  satellites;  instead  it  must  be  modeled.  The  index  of  refraction  is  a  function 
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of  environmental  conditions  in  the  atmosphere,  most  prominently  by  temperature,  pressure,  and 
humidity.  The  total  delay  is  dependent  upon  the  propagation  path  length,  which  in  turn  depends  upon 
the  elevation  angle  ,  of  satellite  j  as  seen  from  site  m  at  time  t,.  The  model  proposed  here  follows 
Goad.'"-* 

The  local  vertical  in  ECEF  coordinates  can  be  computed  from  the  site  longitude  latitude  , 

as  given  by  Equation  (A-28). 

*  cos  P^^sin  +  sin  p^^£  (A-28) 


The  elevation  angle  can  be  computed  from  the  angle  between  the  local  vertical  and  the  site-to- 
satellite  vector  ,  as  shown  in  Equation  (A-29). 

=  -j  -  (r^  •  «^)  (A-29) 


There  are  two  components  in  the  tropospheric  model:  the  dry  and  the  wet.  The  total  delay  is  the  sum 
of  the  two.  There  are  various  parts  to  the  complete  formula  and  these  are  listed  below  for  the  two 
components.  Temperature,  T,  is  in  degrees  Kelvin;  humidity,  H,  is  in  percent;  and  pressure,  P,  is  in 
millibars. 


Dry  Component  (subscript  0) 
Surface  refractivity: 

Height  of  tropopause  (dry): 

WET  Component  (subscript  1) 
Temperature  scale: 

Water  vapor  pressure: 

Surface  refractivity: 

Height  of  tropopause  (wet): 


N„  .  =  77.624 

Omji 


mi 


mi 


0  mi 


5  (0.002277  P  .) 

m.r 

N„  X  10 

Omi 


*Tmi 


7.5  (T^^  -  273.15  ) 
237.3  +  -  273.15  ) 


w  .  = 

X  6.11  X  lo'^"-' 

lif.l 

100 

w  . 

w  . 

N,  .  = 

-12.92  —  + 

371900  — 

lm.i 

T  < 

mi 

tL 

= 

5  (0.002277  )  1 

.  0.05 

J  mi 

1 

O 

X 

T 

\  m.i 

w  . 

m.i 


Expansion  Factors 
Range  to  tropopause 
Where  the  semimajor  axis 


a,  =  6378137  m 
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Factor  a^ 

sin 

'  ~  b.  . 

Factor  b^ 

cosV^, 

ImJ  e 

Expansion  Terms 

Term 

“a 

=  1 

Term  a2 

Term  ttj 

= 

Term 

Term 

=  +  ^bl 

Term 

Term  a. 

Term  a* 

Term 

=  b;„, 

Sum  of  Dry  amd  Wet  Terms 

The  tropospheric  refraction  range  correction  at  each  observation  is  given  by  Equation  (A-30). 


tL  =  10-^ 


(rLf 


(rL/ 


(A-30) 


THE  NAVIGATION  SOLUTION 

At  each  time  of  observation,  a  computed  observation  from  the  receiver  antenna  to  the  satellites  in 
view  must  be  computed.  This  computed  observation,  plus  the  first  term  of  the  Taylor  series 
expansion,  constitutes  a  linear  approximation  to  the  observations,  as  was  shown  in  Equation  (A-1). 
For  single-frequency  pseudorange,  the  computed  observation  (Q  in  Equation  (A-1) )  consists  of  the 
terms  listed  in  Equation  (A-4).  However,  since  two-frequency  data  are  available,  the  ionospheric 
term  can  be  removed  from  the  observations,  as  was  described  in  the  discussion  leading  to 
Equation  (A- 17).  Therefore,  the  computed  observation  has  the  same  form  as  Equation  (A-4)  except 
that  the  ionospheric  term  is  deleted.  The  computed  observation  is  rewritten  as  Equation  (A-31). 


(A-31) 


The  first  term  on  the  right  in  Equation  (A-3 1)  is  the  geometric  range  computed  from  Equation  (A-6) 
and  shown  in  Equation  (A-32). 
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*  iy!i-  (^.i  -  4.)' 


(A-32) 


A  guess  for  the  receiver  antenna  position  (indicated  by  the  tilde)  is  required  in  order  to  compute  this 
range.  The  second  term  on  the  right  is  the  periodic  relativity  contribution,  which  is  computed  as  part 
of  the  satellite  ephemerides.  The  third  term  consists  of  the  two  clock  offsets.  The  satellite  clock 
offset  is  known  from  the  satellite  ephemerides,  while  the  local  offset  is  unknown  (set  to  zero  in 
Equation  (A-31))  and  must  be  included  in  the  state  vector  along  with  corrections  to  the  antenna 
position.  The  fourth  term  is  the  tropospheric  delay  and  is  modeled  as  described  in  a  previous  section. 
Finally,  there  are  the  SA  effects,  which  are  known  in  a  true  PPS  solution,  but  will  remain  unknown 
for  this  demonstration  and  are  ignored. 

Partial  Derivatives 

In  a  navigation  solution,  the  state  vector  (AA^  in  Equation  (A-1))  consists  of  four  elements:  the 
corrections  to  the  antenna  position  Ax„  „  A>'„  „  Az„„  and  the  clock  offset  At„  ,.  The  partial 
derivative  matrix  is  constructed  from  the  partial  derivative  of  each  term  in  Equation  (A-31)  with 
respect  to  each  state  element.  This  matrix  will  be  rectangular,  having  four  columns  (one  for  each 
state)  and  Hq  „ ,  rows,  where  tiq  „ ,  is  the  number  of  observations  (equal  to  the  number  of  satellites 
being  tracked)  from  the  receiver  at  site  m  at  time  t,.  In  order  to  perform  the  navigation  solution,  at 
least  four  satellite  observations  are  needed,  ,  i  4. 


The  partial  of  ,  with  respect  to  x„  j. 

dpt  _ 

X  .  -  x{ 

m.i  -i 

-j 

The  partial  of  ,  with  respect  to^„  ,; 

8|Sl,  . 

-  yj 

sjv, 

'mi 

The  partial  of  p^„  ,  with  respect  to  z„  ,: 

apl,  . 

z  -  z\ 

mi  ‘1 

'mi 

The  partial  of  p-'^  ,  with  respect  to  ct„  ,: 

dpt 

=  1 

d  cx~. 

Equation  (1)  can  be  rewritten  in  matrix  form  with  each  element  explicitly  indicated  as  follows: 
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CD 

SpI, 

siii. 

-J 

^PnLi 

^4..- 

54, 

1 

Pm.i 

P  TtLi 

3p1, 

5p1 

2 

Pm.1 

Pm.1 

^y,.t 

3 

= 

P^ 

r  M.I 

+ 

3p,^/ 

Sffl, 

5p1, 

SpI, 

^4/ 

0 

L  J 

c 

S|il, 

S|Si, 

^4, 

^4/ 

^4* 

Ax 
Ay  . 


Az. 

At 


nti 


m.i 


(A-33) 


Least  Squares  Matrix  Solution 

Returning  to  the  notation  of  Equation  (A-1),  the  solution  for  the  state  vector  proceeds  as  follows. 
Form  the  difference  between  the  observed  ranges  and  the  computed  ranges  as  in  Equation  (A-34). 

—  AX  =  0-C  (A-34) 

dX 


Next  introduce  a  weight  matrix,  W,  which  will  (for  now)  be  an  ,  x  ,  identity  matrix,  where 
Ho  =  j,  the  number  of  satellites.  Then,  pre-multiply  both  sides  by  the  transpose  of  the  matrix  of 
partials  times  the  weight  matrix  as  shown  in  Equation  (A-35).  With  an  identity  weight  matrix,  it  is 
assumed  that  all  observations  have  equal  weight. 

dX)  dx 


AX  = 


dx 


W(0  -  C) 


(A-35) 


The  final  step  is  to  isolate  the  AX  on  the  left  by  pre-multiplying  both  sides  by  the  inverse  of  the  matrix 
quantity  in  brackets  in  Equation  (A-35).  This  result  is  shown  in  Equation  (A-36). 


AX 


V«£ 

-1 

'^1 

dx 

,dXj 

T 

W(0  -  C) 


(A-36) 


The  navigation  solution,  at  the  time  of  the  observations  at  4  is  shown  in  Equation  (A-37).  It  is  the 
sum  of  the  guessed  position  (x^ .,  y„i,  4<)  and  the  deltas  from  the  state  vector  in  Equation  (A-36). 
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The  results  from  Equation  (A-37)  can  be  used  to  compute  an  adjusted  computed  observation  vector 
from  Equation  (A-3 1).  This  adjusted  vector,  p"^,  ,  when  subtracted  from  the  original  observations, 
should  have  a  near-zero  RSS  of  the  mean.  The  standard  deviation  of  the  difference  gives  an  estimate 
of  the  error  in  the  solution. 

In  the  absence  of  SA  (dither).  Equation  (A-37)  would  be  the  equivalent  PPS  solution  that  must  be 
saved.  However,  the  dither  needs  to  be  eliminated  from  the  test  data  set.  This  can  be  approximated 
by  solving  for  the  instantaneous  satellite  clock  offset  in  a  separate  dither  estimator  program. 

Differential  Correction  for  Range  Bias  Estimation 

In  order  to  estimate  the  dither  on  the  satellite  clocks  at  each  observation  time,  at  least  two  static  sites, 
with  known  positions,  observing  simultaneously,  are  needed.  Assuming  independent  solutions,  the 
reason  for  this  is  as  follows.  Each  site  observes  k  satellites,  each  of  which  has  an  instantaneous  range 
bias  due  to  dither.  Thus,  there  are  k  unknowns.  If  the  station  position  and  clock  offset  are  known, 
the  dither  could  be  solved,  but  the  site  clock  is  not  known  well  enough.  It  must  be  included  in  the 
list  of  unknowns,  making  k+1.  A  second  site  observes  the  same  satellites  and  adds  only  its  clock  to 
the  list  of  unknowns.  Therefore,  in  general,  for  I  sites  observing  k  satellites  simultaneously,  there  will 
be  kl  observations  and  k+l  unknowns.  For  k>3,  and  two  sites,  2k  is  greater  than  k+2.  Two  sites 
whose  positions  are  known  will  be  enough  to  estimate  the  k  dither  and  two  local  clock  range  biases. 


The  observation  equations  for  the  dither  estimator  are  shown  in  Equation  (A-3  8).  It  is  assumed  that 
the  site  positions  are  static  and  known.  The  unknowns  will  be  the  two  site  clock  offsets,  and 

the  corrections  to  the  satellite  clock  predictions  found  in  the  broadcast  ephemerides,  x\  +  bx\, 
where  l^j^k. 


P  ^  <  =  r\  +  V  i  +  c  (f  . 

r  m.i  mi  ^  nti  ^  \  ^ 


-  6xf,)  +  Ti, 

-  Ti, 


(A-38) 


Dither  Partials 

The  partials  derivatives  are  as  follows: 
The  partial  of  p^,.  with  respect  to  c  t^.: 

The  partial  of  p^,.  with  respect  to  c 
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Similar  partials  apply  for  site  n  and  the  other  satellites.  The  solution  matrix  is  similar  to  that  in  the 
navigation  solution  described  above.  The  matrix  equation  analogous  to  Equation  (A-33)  is  shown 
in  Equation  (A-39)  below.  The  satellite  dither  solutions  should  be  written  to  a  file  and  saved  for  use 
in  the  PPS  navigation  solution.  The  file  should  include  the  GPS  time  of  week,  the  week  number,  and 
the  satellite  PRN  numbers  with  their  respective  dither  corrections.  The  PPS  program  then  adds  these 
corrections  to  the  clock  predictions  obtained  from  the  broadcast  ephemerides  at  the  particular  time 
of  the  observation.  This  process  approximates  the  dither  removal  that  would  be  performed  in  the 
field  by  other  means  to  obtain  the  PPS  solution. 


As  it  stands,  the  problem  as  formulated  by  Equation  (A-39)  is  ill-conditioned  and  cannot  be  solved. 
In  order  to  fix  this,  one  satellite  was  chosen  to  be  deweighted:  in  effect,  given  a  zero  range  bias.  This 
provides  a  reference  with  which  the  other  ranged  biases  may  be  compared.  Therefore,  the  solutions 
are  satellite  range  biases  with  respect  to  the  reference.  The  desired  range  correction  is  still 
accomplished,  since  the  navigation  solutions  treat  common  biases  as  local  clock  biases.  If  the  true 
local  clock  bias  is  of  no  interest,  this  method  can  be  used  to  remove  the  dither  effect  from  the 
pseudoranges. 
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CONVERT  THE  PPS  SOLUTION  USING  THE  PRECISE  EPHEMERIDES 
With  the  PPS  solution  in  hand,  the  next  task  is  to  improve  it  by  using  the  DMA  precise  ephemerides 
and  clock  estimates  to  replace  the  broadcast  ephemerides  predictions.  To  do  this,  both  ephemerides 
are  required,  as  well  as  the  knowledge  about  which  particular  satellites  were  tracked  at  each  PPS 
navigation  solution. 

FORM  THE  COMPUTED  OBSERVATION  FROM  THE  BROADCAST  EPHEMERIDES 
At  each  PPS  solution  time,  form  the  computed  observation  from  the  PPS  solution  and  the  broadcast 
ephemerides  in  the  same  manner  as  was  described  previously;  see  text  preceeding  Equation  (A-3 1). 
Since  the  data  that  produced  the  PPS  solution  for  this  demonstration  are  known,  use  the  same  set  of 
satellites  that  were  used  at  each  observation  to  obtain  the  PPS  solution.  The  PPS  solution  and  the 
broadcast  ephemerides  will  generate  the  vector  corresponding  to  Cj  in  Equation  (A-3). 

FORM  THE  COMPUTED  OBSERVATION  FROM  THE  PRECISE  EPHEMERIDES 
Evaluate  the  precise  ephemerides  and  satellite  clock  offsets  at  the  same  observation  times  as  the 
broadcast  ephemerides  above.  Use  the  eight-point  Lagrange  interpolation  and  iteration  scheme  to 
find  the  time  of  transmission  for  each  satellite.  Also,  use  the  PPS  solutions  for  the  antenna  site 
positions  when  evaluating  Equation  (A-3 1).  The  result  of  this  will  be  the  vector  corresponding  to 
Cp  in  Equation  (A-3). 

The  relativity  correction  needs  to  be  computed  separately  and  added  to  the  satellite  clock  offset  in 
the  same  manner  as  with  the  broadcast  ephemerides.  However,  since  the  precise  ephemerides  deals 
with  ECEF  positions  and  velocities  rather  than  orbit  elements,  the  computation  of  relativity  cannot 
use  the  same  equation  as  was  used  with  the  broadcast  ephemerides.  An  equivalent  form  is  given  in 
the  equation  below  where  r  and  v  are  the  satellite's  ECEF  position  and  velocity  at  time 


FORM  THE  PARTIAL  DERIVATIVES 

The  partial  derivatives  will  have  the  same  form  as  in  the  PPS  solution  described  previously.  Use  the 
computed  value  from  the  precise  ephemerides  to  evaluate  ,  and  its  components.  At  each  time, 
a  matrix  equation  like  Equation  (A-3  3)  will  be  constructed.  This  is  the  expanded  form  of 
Equation  (A-3). 

LEAST  SQUARES  SOLUTION 

Next,  find  the  solution  to  the  matrix  equation  in  the  same  manner  as  was  described  for  the  PPS 
solution  above.  The  state  vector  LXp  represents  the  difference  in  the  solutions  due  to  the  difference 
between  the  precise  and  the  broadcast  ephemerides.  The  improved  position  is  the  sum  of  the  PPS 
solution  and  the  deltas  from  the  state  vector. 

COMPARISON  TO  THE  TRUTH 

The  improved  position,  as  well  as  the  original  PPS  solution,  must  be  compared  to  the  true  position. 
The  "truth"  will  be  represented  by  on-the-fly  (OTF)  kinematic  solutions  that  are  referenced  to  an 
absolute  point  on  the  main  range  (MBRE).  This  point  is  known  absolutely  to  within  a  meter.  Thus, 
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the  truth  from  OTF  solutions  of  moving  antennas  will  be  known  to  approximately  the  same  accuracy. 
The  error  introduced  by  the  OTF  kinematic  solution  is  much  smaller  than  a  meter. 
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WORLD  GEODETIC  SYSTEDl  1984  (WQS  84} 


STATION 

NAME 


LATITUDE 
(deg  min  sec) 


LONGITUDE 
(deg  min  sec) 


ELLIPSOID 
HEIGHT (m) 


85401 


38  37  39.990N  90  13  19.834W  133.545 


The  position  of  station  85401  is  relative  to  absolute  station 

USAP3  1988 . 


Accuracy  of  the  Positions  and  Ellipsoid  Heights 


The  position  of  the  primary  control  station  USAF3  1988  is 
estimated  to  have  the  absolute  accuracy  of  ±1  meter  (one-sigma) 
in  each  component  (latitude,  longitude,  ellipsoid  height),  with 
respect  to  the  wGS  34  dat'om. 


The  position  of  station  85401  is  estimated  to  have  the  relative 
accuracy  of  ±0.02  meter  (one-sigma)  in  each  comoonent,  with 
respect  co  station  USAF3  1988. 


B-3 


NSWCDD/TR-95/196 


DISTRIBUTION 


Copies 

DOD  ACIWITIES  (CONUS) 

DEFENSE  TECHNICAL  INFORMATION  CTR 
8725  JOHN  J  KINGMAN  ROAD 
SUITE  0944 

FORT  BELVOIRVA  22060-6218  2 

ATTN  CODEE29L 

(TECHNICAL  LIBRARY)  I 

COMMANDING  OFFICER 
CSSDDNSWC 
6703  W  HIGHWAY  98 
PANAMA  CITY  FL  32407-7001 

ATTN  B  ROTH  I 

SMALYS  I 

DEFENSE  MAPPING  AGENCY 
4600  SANGAMORE  ROAD 
BETHESDA  MD  208 1 6-5003 


NON-DOD  ACIWITIES  (CONUS) 

ATTN  GffT  AND  EXCHANGE  DIV  4 

LIBRARY  OF  CONGRESS 
WASHINGTON  DC  20540 

THE  CNA  CORPORATION 
PO  BOX  16268 

ALEXANDRIA  VA  22302-0268  1 


INTERNAL 


E231  3 

E282  SWANSBURG  1 

K  1 

KIO  1 

K12  1 

K12  HERMANN  5 

K12  SWIFT  1 

K13  1 

K14  RISINGER  5 


(1) 


